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ABSTRACT 

A finite-difference solution for laminar viscous flow 
through a short wavelength sinusoidal symmetric as well as non- 
symmetric channel and through a converging-diverging pipe is 
presented here^ The physical wavy domain is transformed into a 
rectangular computational domain in order to simplify the 
application of boundary conditions on the walls# The discretized 
conservation equations for mass* momentum and energy are 
derived on a control volume basis# The resulting discretized 
equations are solved using the SIM.PLEG method suggested by 
Van Doormaal and Raithby [12]. We also follow their other 
suggestions regarding i) the use of E-factor in place of the 
usual under-relaxation factor^ and ii) the solution of the 
pressure correction equation. The latter was found to be very 
effective in reducing the number of iterations as compared to 
the procedure suggested by Patanlcar[ 7]. Results are obtained 
for both the developing and the fully-developed flow for a 
Prandtl number (Pr) of 0 * 12 , for the Reynolds nuraber (Re) 
ranging from 100 to 500, and for non-dimensional wall arrplitude 
(X/l) = 0.10 to 0.25. Solutions are obtained for constant 
fluid properties and for a prescribed wall enthalpy only. It 
is found that the separated flow region grows with increase 
of Re and X/L and is mainly confined to the diverging part 
of the symmetric channel or pipe. However* for higher Re and 
X/L it occurs in the converging portion as well. In the non- 
symmetric channel flow separates near the upper boundary in one 



x±x 


section whereas it s^arates near the lower ■ boundary in the 
adjacent section* Per-cycle pressure drop in the fully- 
developed region is same and increases with increase of Re and 
X/L* A point of inflection in the enthalpy profiles is observed 
in the s^arated region whereas in the non-separated region the 
enthalpy distribution is almost parabolic* 



Chapter 1 


INTRODUCTION 

A lot of effort has been expended during the last few 
years to enploy boundary fitted coordinates in conjunction with 
finite-difference techniques to solve fluid flow problams in 
and around wavy surfaces- The analysis of such problan finds 
application in different areas such as generation of wind 
waves on water/ formation of sedimentary ripples in river 
channels and dunes in desseibs/ transpiration cooling of re- 
entry vehicles and rocket boosters, stability of liquid film 
in contact with a gas stream, rippling of melting surfaces, 
cross hatching on ablative surfaces and film vapourization in 
combusiiovi chambers* This type of configuration is also used 
in heat exchangers in order to enhance the convective heat 
transfer characteristics* Physiologists are also interested 
in it trying to explain blood and urinary flow and to apply 
the results for optimal design of artificial organs* 

Earlier development 

One of the earliest studies on viscous flow in 
sinusoidally varying channels and pipes is that of Bums and 
Parks [1]* They carried out the solution for small amplitude 
channels and pipes by expressing the stream function in a 
Fourier cosine series under the assumption that the Reynolds number 
is small enough for the Stokes approximation to be valid* 
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T sangarls and Leiter E 2]| solved the same problem by expressing 
the stream function in a Fourier series not in the physical 
plane but in the transformed plane where the wavy boundary is 
transfomed into a straight one* Their above proposed analytical 
perturbation method for creeping flow was extended to laminar 
flow for higher Reynolds numbers by the same authors [ 3] and 
results are presented for an amplitude to wavelength ratio of 
0*1/71 • Steady and unsteady flow through furrowed channels has 
been studied numerically by Sobey [ 4 ]. He especially 
concentrated on the Reynolds number effect for separated flow* 
Fluid flow connected with heat transfer in wavy channels was 
calculated by Vajravelu [ 5] by the perturbation method for long- 
wavelength channels* Heat transfer effect has also been studied 
by Prata and Sparrow [6] for laminar flow in an annular duct 
having a streamwise— periodic variation of the cross-sectional 
area* They obtained the solution by using the SIMPLE algorithm 
of Patankar[7] for Reynolds number ranging from 50 to 1000 and 
Prandtl number ranging from 2 to 10* The viscous flow in a 
circular pipe whose radius has slight sinusoidal variation was 
treated by Bel infante [s] using the power series method# Some 
experimental work has also been carried out by Hsu and Kennedy [9] 
to find the variation of pressure and shear stress along a wavy 
pipe for turbulent non-separated flow. Stephanoff^ Sobey and 
Bellhouse[lO] Compared the results of Sobey* s calculation with 
e 35 >erlmental observations by visualization techniques. 
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Present Work : 

In the present worlc, a finite-difference solution for 
two-dimensional, viscous laminar flow in sinusoidally varying 
channels and pipe is presented* To solve such problOTs, both 
orthogonal and non-orthogonal grids have been generated numerically 
in the physical domain so that the boundaries are grid lines* 

The governing transport equations are then solved in the 
generalized coordinate system* In the present development, 
however, suitable non— orthogonal algebraic transformation are 
used in order to transform the physical wavy domain into a 
rectangular computational domain* The discretized conservation 
esquations are then derived on a control volume basis [ll]^ This 
approach has two attractive features* One is that it facilitates 
physical interpretation of the terms that result from the 
coordinate transformation* Thus for example pseudodiffusion 
terms that result from the non-orthogonal nature of transformation 
can be easily Identified. The second attractive feature of the 
control volume is that it ensures global conservation of mass, 
momentum and energy* 

The resulting discretized equations are solved using 
the SIMPLEC method recently developed by Van Doormaal and 
Raithby This method is superior to the SIMPLE or SIMPLER 

method suggested by Patankar C? j. It is consistent interms 
o f approximation made and avoids the task of searching for 
optimal values of under relaxation factor since in this method 
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the relaxation factor for pressure is simply unity* Their other 
suggestions regarding use of E-factor in place of the usual 
under-relaxation factor and the solution of pressure correction 
equation are also adopted here* The latter was found to be very 
effective in reducing the number of iterations required for 
convergence as compared to the procedure suggested by Patankar* 
This reduces the solution cost considerably* 



Chapter 2 

ANALYTICAL FORMULATION 


Description of the Problan : 


We are concerned with two-dimensional/ viscouS/ laminar 
flow through a wavy channel and a wavy pipe formed by altematinc 
converging and diverging cross-sections defined by a sinusoidal 
function as well as through a wavy channel formed by a 
sinusoidal function but having a constant cross-sectional 
area in the flow direction* For all the above three cases, the 
solid boundaries are curved and do not lie along a coordinate 
line* A schematic view of the types of physical domains being 
considered is shown in Fig* 1* The walls of the channels or 
pipe are denoted by a sinusoidal function* For the non-symmetric 
channel the sinusoidal function is different from that for the 
symmetric channel or pipe# These are given by 


For the symmetric channel or pipe 

6 "(Xj) = L - X(1 + cos (71 X^/L)) (2.1*1) 

For the non-symmetric channel 

6-'(X^) = X Sin (ti X^/L) (2*1.2) 

whore L is a reference length and X is the amplitude of 

the channels or pipe* The origin of the coordinate system 
iat placed on the entrance plane of the duct as shown in Fig^ 1, 
Here X^ and X^ are used as the generalized coordinates* For the 
Cartesian coordinate system used for channels Xj^ and X^rtpresent tl- 





L+6(Xi) 



X2= 6(Xi) 


(a) Converging and diverging 
symmetric channel 

(b) Non symmetric curved 
channel 


j — X 2 = 6 '(Xi) 


K) 


(c) Converging diverging pipe 


Fig. 1 Types of channels and pipe considered 
(one cycle only) 
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X- and y~ coordinates respectively, whereas in the cylindrical 
coordinate system used for pipe represents the z -coordinate 
and X 2 represents the r-coordinate» 

Other typos of problans, like external flow over a wavy 
boundary can be dealt with by the present methodology tut these 
will not be considered here* 


2*2 Conservation Ecruations : 


The governing eguations to be considered here are the 
continuity, momentTjjm and energy equations# Constant thermo- 
physical properties are assumed and viscous dissipation is 
neglected in the energy equation. Moreover, the flow is either 
two-dimensional or axisymmetric, i.e#, the flow parameters 
are independent of the z— direction in the case of a channel and 
independent of the ©-direction in the case of a pipe# Using 
these assumptions the governing equations become : 


a) 


For channel : 


Continuity equation : 


3U. 


X. 


su 


X, 


3X. 


+ 


ax. 


o 


( 2# 2#l) 


X-i -momentum equation : 

•• -r-r TT 


au 


X. 


au 


X. 


u 


X. 


ax^ 


+ u. 


1 ap 


+^-. ( 

p ax. ^ P ^ 


X 2 ax^ 




X. 


a V 


X, 


ax: 


ax^ 


( 2 # 2 # 2 ) 


X 2 *-momentum equation 


au. 


dU 


X, 


u 


Xl 9^1 


+ u. 


Xg 3X2 


1 ^ M 
p 9X2 P 


1^, 


X. 




X. 


0 ( 2 . 2 . 3 ) 


ax‘ 


ax: 


2 
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energy equation s 

% ^ 3 X 3 

b) For pipe : 


JL M 

PC '■ 2 2 '' 

^ P 3X^ 3 X 2 


(2.2*4) 


Continuity equation : 
3U 3 U 

^0 ^0 
-mom entum equ at ion : 


( 2.2.5) 


3U 


X. 


3U 


X. 


u 


Xl 3^1 


+ u 


X 2 3 X 2 


i _ - i 9£_ + tL (. 


3^^ 3 ^, 


X, 


P 3 X^ ' P '9^2 


3“U 3U 

1 + i + -1 -i) 


3^2 


X 2 3 X 2 

( 2 . 2 * 6 ) 


X 2 -mom entum equation : 


3U 


X, 


3U 


X, 


aV 


U 


Xl 9X, 


+ u 


X 29 X 2 


1 !£_ + ( 

p 3X2 P 3 x 2 


X, 


a^ 


X, 


3x; 


X 2 3 X 2 


^2 -Xg 

-^2 


) 


(2*2.7) 


energy equation : 


-U^ IS- + U 


3^ . 3^h . 1 3h’ 


■p^” 2 ^ 2 

3XJ 3 X 2 ^^2 


X^ 3X^ ■" -X2 3X2 " PC^ \^2 ”^2 x‘, 3X,^ 


( 2 . 2 * 8 ) 


Boundary cwnditions for the above problexis will bo discussed lateri 


2.3 Formulation For Symmetric Channel and Pipe t 

Since both the converging-diverging channel and pipe 
are synmetric about their respective centre lines, the total 
solution domain can be folded at the centre line with certain 
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modifications in the boundary conditions to get the numerical 
solution of the problem* The coordinate transformation described 
later for this channel and pipe is also of the same type* 

This is made possible by introducing a factor K which is unity 
for the Cartesian coordinate system and ^two' for the cylindrical 
coordinate system in such a way that the same mathematical 
analysis can be applied to both the cases* 

Therefore with the insertion of factor K, the governing 
equations take the form 


Continuity equation : 
3U 3U^ U„ 

A- A 

dx^ 3 X 2 X 2 

X^-momentum equation : 

a,u. 


( 2*3*1) 


u 


Xi 


4. u - 1 + a (V^J ) 

^X^ 3X^ P 3X^ P ^ ^X^^ 


( 2*3*2) 


‘2 ^ 2 

X 2 -momentum equation ; 


8U^ 3U^ 


U 


u 


X 


1 


+ u 


X 2 3 X 2 


2 _ 1 3P U , ^ 


X, 


X, 


energy equation : 

2 

where the La^lacian (v ) operator is 
v2 ^ + ( K-1) |r^ 

3x“ ax| ^2 ^ 2 


(2*3*3) 


( 2.3*4) 


( 2 - 3 . 5 ) 
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The boundary conditions (except those at exit) are 


( 0 ^X 2 ) = (assumed constant here, although a 


function of X 2 


ax, 

Ux^(Xj, 6 VX^)) = 0 , 

Ux (Xj_,0) = 0 

Ux^(Xj,5 

p(0,X2) = Pq 

h(0,X2) 

h(Xj,6'(Xj)) = h^ , 

2.iL. ( X 1 

8 X 2 ^ 1 ' 


( 2-3.6 a) 


where and h^ are the velocity, pressure and enthalpy at 

entrance to the duct and h^ is the enthalpy at the wall. 

Boundary conditions at exit 


The characteristics of velocity and temperature (or 
enthalpy) in a periodic fully developed regime are discussed in 
D6]- The periodicity conditions of velocity can be expressed 
as 

(2.3.6b) 


U^^(x£,X2) =U^^(xJ~VX2) 


where 5^ is an arbitrary station in the fully developed region 
and is the cycle length- 

For a periodic thermally developed regime with constant 
wall enthalpy, we have 


H(3^,X2) = H(3q-L^,X2)‘ 


( 2#3.6c) 
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if the dimensionless enthalpy H is defined as- 


h - h. 


H = 


w 


hj^-h 


w 


wh 


ere the balk enthalpy h^^ is given by 


»’b = »'w + 


A 1 

/ U„ dA 
A 1 


and 


(b) 




H ( Xj X2 ) H ( X2) 


( 2 • 3*6d) 


if H is defined as 


h - h 


H = 


w 


h^ “ h 
o w 


Now using the following non-dimensional variables 

u. 


X X '^x 

X Y = -^ , u 


X, 


, V 


u 


h - h. 


H 


w 


^o “ ^ 


, p 


p - p. 


o 


pu‘ 


the governing equations take the form 


( 2*3*7) 


S + |y + Y = ° 


V 


TT O. ^7 


u 


dV 

dx 


dv 


( 2 * 3*8 } 


+ V ^ s 


ax dY ^ 

- i + A ^ ^ i Iy- 

( 2 * 3 . 10 ) 
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( 2» 3 mil) 


PU L 

where Reynolds number/ Re = — —— and Peclet number, Pe = Re#Pr; 
Pr =MCp/k. 

The boundary conditions become 


U(0,Y) = 1.0 
|y (X/0) = 0.0 


U(X,6(X)) = 0.0 
V(X/0) =0.0 


U(X*,Y) = UCX'**' - L*/Y) 

C 


V(X/6(X) ) = 0.0 / 

V(X’'^,Y) = V(X* - lJ/Y), P(0/Y) =0.0 
H(0,Y) = 1.0 / 

|S(x,o) = 0.0 , .X) = 

^ ^ H ( X'' Y) H ( X*-2L*/ Y) 


( 2.3«12) 


H(X, 6(X)) = 0.0 


H(X*-L*/ Y) 


X 


.* 


6'(X) 


where x'^ = ^ , lJ = L^L and 6(x) =2 

Equations( 2.3.9 ) and (2.3.10) are the X- and Y- components of 
the vector momentum equation. Therefore, the momentum equations 
(2.3.9) and (2.3.10) may be expressed as 


lx + + 


^ - 1 (iiu . ^ , /K-1) I e 

lx Re ^ ^ + (K 1) Y ej^ 


+ {U + V -^Y + ^ Re + ^„2 + ^ 


ax 


aY‘ 


(K-l) -4)} = 0 

Y^ 


( 2 . 3 . 13 ) 


where e^^ and e^- are the unit vectors in the X- and Y~directions 
respectively. The need for writing the momentxim equation in 
vector form will be justified later. 
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The next step in the analysis is to introduce a 
transformation of coordinates in order to transform the non- 
rectangular physical domain into a rectangular computational 
domain so as to avoid the task of numerically generating the 
boundary fitted coordinates and also to simplify the application 
of boundary conditions at the walls- 

Specifically,^ the X and Y coordinates are transformed to 
V and g coordinates respectively by the relations 

= X , g = Y/6(X) (2.3.14) 

such that I = 0.5 at all po'ints on the curved boundaries. 

The lines of constant T? and § for a given boundary 
shape are illustrated in Pig. 2a. It is evident that a control 
volTome contained between lines r) = Tjy “ ^2 ^ “ ^1^ ^ ~ ^2 

is a curvilinear element with non -orthogonal sides* The 
quantities e ^ and are unit vectors in the physical cooirdinate 
system, which lie along the lines of constant g and constant t? 
respectively. It is evident from the figure that the direction 
of e^ changes with position. This is in contrast to the unit 
vectors e^ and which do not change direction throughout 
the solution domain. The shape of the control volume element 
both in Cartesian coordinates and in cylindrical coordinates 
is shown in Pig* 2c and Pig. 2d respectively- Now the vector 
momcaitum equation (2* 3*13) is recast in terms of the unit vectors 
e^ and e^ and the coefficients of the respective unit vectors 
are identified as the momentum equations in the 1) and | 



(a) Shape of control volume 
in physical space 



(c) Shape of control volume 
element In cartesian 
coordinates (dz = unity) 


(b) Shape of control volume 
in h “5 plane 



(d) Shape of control volume 
element in cylindrical 
coordinates (d0 = unity) 


Fig. 2 Shape of control volume for a symmetric 
channel or pipe . 
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direction^. These momentum equations are used to r^lace the 
X~ and Y-inomentum equations- The continuity and energy equations 
are scalars and remain as they are* 


In order to resolve the momentum equation along ^ , g 
directions, the unit vectors e^ and e^ illustrated in Pig* 2 are 
used* Since lines of constant t) coincide with lines parallel 
to the Y~coordinate, we have 


eg — &Y (2*3*15) 

To determine e^, consideration is first given to the unit 
vector n which is normal to a line (or surface) of constant § « 
Since the gradient of g is normal to a line of constant g , 
we have 

^ grad^ vg 

^ 

IgradS I ivgj 

-* 9g -♦ 3 g 

Now Vg = e^ + Gy- ^ • 

From the definition of g (Eq* (2*3*14)), we get 


and 


M _ 1 

3Y " 6 

Yd6_ id6_ 8 

■53^“‘"g2dX"~6^n""’6 


where 


0 



vs = I (-P 


( 2*3*16) 


7*8 e^^ + Gy *"8 ®x ®Y 
'/■(It 8^) 


and 


(2*3*17) 
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where a = 1 + |3^ (2.3-18) 

Since the unit vector e^ is perpendicular to the vector we 
have 


e^n = 0 ^ which gives 
_ ®X ^ ®Y 

s -p = j^y' 2 f 2 • 3 • 19 ) 

This unit vector is not constant/ but it is a function of positioni 

The following inverse relations can be obtained directly 
from Eqs» (2-3*15) and (2-3-19). 


ey = 

Now any vector a can be expressed as 


a = e„+ a^ or a = 


SO that 


ay = a.Oy = (a e t a e ) -"^ 


V V £ S 

[a„ (■ 

/o 

13 -h 


^ \ ©jr 4 “ 3 1 * 


and 


®X “■ 


a-Ojr = (a^ e^ + a^ e^.e^ 
r = [ a ^/^( ^ 


a^a 


- 1/2 


( 2-3.20) 


(2- 3. 21a) 


(2-3-2lb) 
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From Eqs* (2*3»2l) we get 


= a 

Tl 


1/2 


- "3 

Therefore/ ifa=U = U &^+V^=U^e^+U^e^ 

X x ri n ^ ^ 

Then = a^/^' U 


( 2.3*22) 


Ufc = V - 3U 


( 2*3-23) 


Substituting and eyin terms of e,^ and e^ from equation 
(2*3-20), equation (2*3*13) becomes 


3U , „ 3U 3P 1 ,3% , 3% 

9 9 


+ ( K ~ i ) Y 3 Qg ) 


4 - fTT 4. V 3V . 3P 1 /9^v 3^V ^ 1 3v Vs - ^ 

+ ^ ^ a 7 + -^Y “ Ri + -^ + (K- 1 ) Y gY- (K-l)^)}e^ =0 


or 


rniU - v ^ + ^ P - A . 4. 9^ . ( t , .s 1 aUs j 

{U 35 ^ + V p. + Re ^“2 + 771 + i; Y ^ *=77 

9 9 X 


,1/2 




ft M.) 4- vf ^ - R 2^^ - fft ££ - £ i ;7- 4. 

3 3 gY ^ ^3 gY Re \^2 ^ 


3 V 


3 U , 


3P 3P^ 1 /3^v . a^V 


ax '* 3 Y 

+ ( K -1) 1 1^)1 3, = 0. 


2r. g2 

Y*' 3X'“ 3Y 

From this vector equation it is possible to write 


i) momentum equation in the 7>-direction ; 


3 U 


U^+v|^ = -|y + -4 ■*■ 7^ ■*■ i (2*3-24) 

d X 3 i 

ii) momentum equation in the f -direction : 

,,/3V „ 3Us ^ ,,,3V . 3Us 3P 3Ps ^ 1 

U (^ “ 3 ^) + V(p - ^ p ) = (0 p - p ) + ^ 


3P 


1 rd'^V 




1 3 Us 


(g - g - CK - a ) i I ? - ( K - a ) ^ ( K - i )| ||), 


( 2-3*25) 
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Since V = Ug + |3U, we get 

aV , rr dB ^ 

alx=-^ + Ug.|+3-3^, 

^ + IT M 4. o 


h. B 38 3U 

- ax ^ 0 

_ , au 90 ^ tt 

- 2 ^ • si + ^ 


3X 


2 + ^73 


a^ ^ 

*T" 


% 


and 


B^'^ ^ jy. 30 . TT 3^0 . « a 


3^1 = 2 1? • If + u ^ + e 


3X ax“ 

h 


2 ■t o 


Using these derivatives/ ^ -moment'Um equation becomes 

If + If + ^) = 3 H - if 


1 r ^ SU Be , rr B^B . ^ B^ 




9 

^ .apJa.. 


^ + 2 ^ ^ 


ax‘ 

a^Ut 


ax 


ax 


n ^ *1 *s n <\TT 9^fc 


aY 


or 


*• (k— 1 ) (u^ + 

b 

0U) 

aut 

+ V 

II 

^|n 

<0 fco 

^ H 

- U(u II 4 

■ '^If’ 

A 

+ 2 gr3? §T 

■*■ 2 |y 

|| »< 

The form 

of Eq» 

(2.3 


■ . a^ 

7" ' ' ' A *T“ '■"' '■^""■"yt 

ax aY"^ 


1 au. 


a 


ax' 


a^t 


ax^ 


aug^ 

^aY - 


ax ax-' 


(2*3*26) 


U 


g as the primary dependent variable instead of the cltunsy 
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corabination of U and V that appears in the convective and 
diffusive terms of Rj. ( 2.3*25) ♦ On the other hand Eq. (2.3*24) 
wherein U plays the role of the dependent variable/ is much 
simpler than the transformed version having as the dependent 
variable- Therefore it is advantageous to retain Eg, (2-3*24) 
and to use U rather than U^- 

The boundary conditions for Eqs- (2-3-8)/ (2-3.24)/ 

(2-3.26) and (2-3-11) are 

U(0/g) = 1,0 / U(n,0.5) = 0.0 / H (71/0) = 0.0 

U(1?^^) = U{( 71*~ L*)/ ^ U^(0,g) = - 3^ U(0/|) = - 3^ 

Ug(71/0) =0.0 / Ug(T)/0.5) =0.0/ U (71 

P(0/g) = 0.0 / H(0/g) = 1-0 / H(7), 0.5) = 0*0 , || (TJ/O) = 0*0 / 

H(T)*g) H{(71 *-l*) /^} 

H{ ( 77*-L* ) , g } H{ ( 7]*-2Lp / t } 

where 77 ^= X* and 3 ^ is the value of 3 at entrance to the 
channel or pipe. 

Control volxome form of the conservation ecaiations i 

At this point/, the governing equations [Eqs* (2*3.24), 
(2-3.26), (2.3*11), and (2.3*8)] are integrated over a control 
volume in physical space bounded by lines of constant 77 and 
constant § • Such a control volume is Illustrated in Fig, 2a. 

This is done with a view to obtain the discretization equations. 
Equations obtained in this manner express the conservation 
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principle for the dependent variable for the finite control 
volume, just as the differential equation expresses it for an 
infinitesimal control volume* Thus the resulting solution 
will exactly satisfy the integral conservation of quantities such 
as mass, momentum, and energy over any group of control volumes, 
and thus over the whole domain* 


The term-by-term integration of Eq* (2»3*26), using the 
divergence theoranji, will now follow* 


3U. 


3U, 


; (u + V d¥ = ; (u.vug) dv = / (u.n)Ug ds 


( 2,3*27) 


and 


, aUt 

/ ( ^ 4. - 4. (K-l) « — Y ■ ) d¥ 


V ax' 


3Y' 


S d¥ = ; V*VUf. d¥ = / n.vu. dS 

V ^ V ^ s ^ 


(2.3*28) 


With these, Eq* (2*3*26) becomes 


/ (iJ.n)Ug dS - ^ (n*VUg) dS = / H 

- / U(i5-.V3)d¥ + ~- / {Uv% + - (K-l) (Ut+j3U)} d¥ 

y- K 0 yr y S 

( 2*3*29) 

Using the vector identities similar to those in Eqs* (2*3*27) 
and (2*3*28), Eqs. (2*3*24), (2.3*11), and (2*3*8) become 

/ (t)*n)U dS - ^ ; (n.vU) dS = - / ^ d¥ (2*3*30) 

S S V 

/ (U.n)H dS - J (n.vH) dS = 0 

S ^ s 


(2.3*31) 
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I (n.U) dS = 0 (2-3*32) 

S 

where V and S represent, respectively^ the dimensionless volume 
and surface of the control volume* 


For evaluation of the surface integrals ej<presslons are 
needed for the area dS of the surface^ the gradient operator 
and the unit vector n* To derive these quantities# it is 
necessary to consider a formal transformation between derivatives 
with respect to Y and those with respect to T) # g where r) and 
t have been defined in Eq- (2-3*14)* This transformation is 


3 3 dV d 3f3 B3 


(2*3*33) 


3 3 S'H 3 3l 1 3 


(2*3*34) 

The surface integrals v\^ich appear in Eqs- (2 * 3 * 29) — (2 - 3 *32) 
arc evaluated with reference to the control volume shown in Pig* 2< 
Also the total surface intesgral is subdivided into a sum of four 
surface integrals over the segments S 2 # and S^* 

The first surface integral in Eq. (2*3*29) will now be 
ovaluatod- 


For surface 

Unit normal vector n = e^ 

Therefore# U.n = (U e^ + V = U (2*3*35) 

Elemental surface area dS = Y^ ^ dY* 
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Now from Eq* (2*3*14), Y = 


m» dY = ( since 77 is constant along S^) 


Thus dS = ( 5 g) 


For surface S, 


K -1 


6 d| 


{ 2*3*36) 

(2*3*37) 


^ *“ 1 / ? /- ^ 

Unit normal vector n = cc + eyj 


U *n = (U + V ©Y") • ( "“3 t ©y) ^ 


- 1/2 


,- 1/2 _ 


= Ufc a 


= (-3 u + v)a' 

Elemental surface area 

dS = Y^“^ (dX^ + dY^)^/^ 


- 1/2 


(2.3*33) 


Since rj = X and Y = | *6( T]) , therefore dT? = dX and 


^ • d77 ( V I is constant along S^) 


3 dn 


Thus dS = (ag)^"^ (1+3^)^^^ drj= (6g)^**^ (2*3*39) 

For surfaces and dS is identical to those for and S^/ 
with the exception that the outward normal n has the opposite 


sign) 


Thus, the first integral in Eq* (2.3*29) is 


/ (t).n)UtdS = / (0*n)Ufc dS + / (U.n)U dS + / (U*S)U dS 
S ^ ^ ^2 ^ ^3 ^ 

+ / (U.n)Ug dS = / UUg(6g)^'"^ 6 <ag 


s 

i 

+ S UfeUfcXaS)^*”^ dT7 - / UU (61)^"'-^ 6 dg 

c s S S s 


,K -1 


-/ UtUfc(6g)^““^ d’) 

S4 


(2.3*40) 
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From the above. It is observed that the term S-n yields only 
one velocity component* This is because the unit vectors 
and eg were chosen to be tangent to the coordinate lines and 
the velocities were resolved in those directionsf otherwise x 5 #n 
for surfaces S2 and would have involved both velocity 
components U and V. 

The V operator appearing in the integral has its 
standard form in X-Y cooirdinates as 


V 




- 3 
®Y aY 


(2.3*41) 


3 9 

which, after substitution for and from equation ( 2 * 3 * 33 ) 
and (2*3*34), becomes 


V = 


3 ^ 6 ar 


(2*3*42) 


^ Q 


Therefore, for surface S^, 

n*Wg = * (^ ®x + "p" ®Y^ 


3 X 


317 


6 35 


For surface S 


2 ' 


3 U, 


3 U 


n.vu 


= CX ^/^( <“8 ©j^ + e^) gY ' ®y^ ” ^ 


'3- 

„-v 2 , ^^5 * 1 !!ii) . 

“ 3fr * s jf * z-sr’"~ 


( 2 » 3 ^ 43 ) 




3 X 




3 § 


(V a » 1 + 32 ) 


(2*3-44) 
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Aside from a minus sign^, n.vUg for surfaces and is 
identical to those for and S 2 respectively* 

With the aid of the above, the second surface integral of 
Eq# (2* 3*29) may now be expressed as 


/ (n.vU ) dS = / (n.vu-j.)dS + / (n.Wj.) dS + / (n*vu.) dS 
S ^ S, S S, ^ ^ 

123 

+ / (n.vu^) ds = / eas+/ 

4 1 '^2 

- 3 gif) (6 ^) “ dr? ~ / (— - i -^){ 5 g)K-l 5 

3 

,CL 3Ufc X—1 

” s 5 ~di ” ^ 3^ ^ (2*3»45) 

^ aug 3Ug 

In the above equation the terms 0(g^) and are cross- 

diffusion terms across surfaces or and S 2 or driven 

respectively by velocity gradients in the § and r? directions* 

They originate from the non-orthogonal nature of control volume 

boundaries that appear through the transformation given by 

Eq* ( 2 • 3»33) • 


Attention will now be focussed on the volume integrals 
that appear on the right side of Eq* (2*3*29)* 

The volume element d¥ may be expressed as 

d¥ = dY dX = 6 d§ d r? (2*3*46) 

Now various terms related to the right side of Eq* (2* 3*29) 
are evaluated as follows i 
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Since 3 = 


and = 


Thus U.vp 


= (U ej. + V ■^) ,(|| + = ^|l+^a^* 

^ = 6 in ' therefore || = i 


Y d<5 do . fc d^6 _ ^ d^6 _ ^ 


_ y CIO do d 6 

p 35 i • OT + i ^ = s 


dn 


2 6 0 ?) 


= u€-^^- 2 JM + v d 6 _i da , d ^6 

^ ^ 7,2 6 dl? + 6 3 n - 6 d^ ^ 


u 




( 2«3-47) 




Now ^ 
3 X^ 


J-| + + (K-1) I 34 

3 X 3 Y 3 

^ ^ - 1> f - 'h ' '-!!># 


dr) 

,3 


dr)" 


■ d »)3 62 W 6 ' 


Similarly 

9^5 


• 2 _ 

• • V 3 


|y (|f) =° 


a avi' 


_i_ 1 da i_ da 

a ^ " £^2 a?) 


fc 0^5 ^ 2 | /das 3 


^ 4f^+ (K-l) (2*3,48) 

^ df) ga^ 


VU,V 0 


= ^ + |t *^§1 ®X + fl 


ev) 


3^ 3B j. 9U 3 

3 lt -ft + ^ *3 
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fdu 
" ^3r) 

! It) <5 

sag ^ 

^ 6 - 
^ 2 
d’T) 

£ J 
6 « 

^) + 1 2 H 

^ 6 3g 

1 ^ 

' 6 S) 



= ^ (1 
- a 7? 

dT?^ ^ 

~) + 
dT)^ ^ 

^ (l±af ^ 

a§ ^52 d77 

- M. afaj 

<5 



= 2 .u (. 

- an ^ s 

d77^ 

^) + 
dV^ + 

ag 

f d 6 _ , 

\2 dv 

0 

ei £6) 

^ dri^ 

( 2 , 

» 3 * 49 ) 

^ dX 


- i 

5 ag^ 

- 

1 _ fi 

S 3§ - ^ 

^ - a ap 
ar? 6 al 

( 2 . 

. 3 * 50 ) 

and 








-ii_ (u 

y2 

+ I3U ) = 

sV 

+ 3u) 


( 2 . 

.3*51) 


All the ingradients needed to evaluate the right side of 
Eq. (2*3»29) have been derived* Therefore^ their sum/ denoted 
by b/ is giv<a:i by 


b 


U ' 2 

- ^ m + ^ 6 d§ dn 

V 5 dn^ 

aaSdii+X;^(|u ^ - f 

+ "Ir ^ ^ (65)^“^ a a| dv 

3 5 OT) 6 caT)^ 

- { (K-1) (Ug + PU)} (6t 6 dg dn 




(2.3*52) 


Now with the aid of expressions in Eqs* (2*3*40)/ (2*3*45)/ 
and (2*3.52), Eq#( 2*3*29) becanes 
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. aUfc 3 Ut 

-j fuuj 5 -i{a^-P 5 ^)! (61)="-^ as 

3 

■1 a K-1 

+ f^S^S - ^ (f -9! - p ^ 

2 

- „ 3Ufc aufc ^ . 

- i ^ 'f -iF - P -sTr'i ( 3 77 = b (2.3-53) 

4 

Because of similarity among Eqs* (2*3#29) - (2.3#3l), it is 
possible to rewrite them in a common compact form by introducing 
the following parameters : 


Q = - 


a 

6 


dJP 

ag 


A= a 


dV 


y = 


dV 


5? = 3 


2JP 

di 


( 2»3.54) 


where <p stands for U,Ug or H< Using the above parameters^ the 
common form of Eqs. (2-3*29) - (2.3*31) is 


/ £U<P6 +P(y+f)| (6g)^”^ dg - ; {U(P6 +P(r+^} (6§)^~^ dg 

1 3 

+ / {U|<P +r(Q+A)} (6 g)^-^ dT? - ; £Ug«P + p(Q+A)} (6g)^”^ dT7 

^2 ^4 

= b (2*3*55) 


where P is the diffusion co-efficient and is equal to 1/Re for 
momentum equations and 1/Pe for the energy equation* The b- 
term for g -momentxam equation is given by Eq* (2*3*52)* For 
7) -momentum equation, b is obtained from the right side of 
Eq. (2-3*30) and is given by 
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b = -/ ||d¥=-; 6 dV (2*3.56) 

y: ^ dV 6 5^ 

For energy equation [Eq. (2*3»3l)] b is simply zero* 

The terms U and in Eq. (2*3*55) are usual 

convective fluxes crossing the control volume surface* The 
quantities 7 and q are, respectively^ proportional to 
diffusive heat or momentirm fluxes in the TJ and | directions* 
The i\ and ^ terms can also be thought of as representing 
pseudodiffusion fluxes* ThuSj, A denotes diffusion driven 
across a § -f ace by a gradient in the ^7-^irection, while W 
represents diffusion across an Tf-face driven by a | -gradient* 

Finally, the continuity Eq* (2*3*32) can be written as 


J 

S. 


(n.t?) dS + / (n*u)dS + / (n*u)dS + / (n.if)dS = 0 

S„ S. 


or 


I U(6|)^"^6 dg + / Ufc( 6g)^'^ d?l - / U( 6g)^“^ 6 dg 
"l "2 

- / Ufc(dl)^“^ = 0 (2*3*57) 

^4 

If velocities are retained interms of physical components U 
and V, the resulting continuity equation would have involved 
two extra mass flux terms on surfaces S 2 and S^* 

2*4 Formulation for Non-svtnmetrlc Channel 


Using the same assumptions and same non-dimensional 
parameters the governing equations are the same as earlier 
f Eqs * ( 2 * 3 * 8 )# (2»3*9)# ( 2 * 3 * 10 ) and (2*3*11) with K =s ij* 
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The boundary conditions in this case are/ 

U(0/Y) = 1.0 / V(0/Y) = 0.0 

U(X^ 6 (X)) =0.0 / V(X, 6 (X)) =0.0 

U{X,(1+6(X))} = 0-0 / V{ X/ (1+6(X))} =0.0 

(2.4.1) 

H(0/Y) =1.0 
H(X^ 6 (X)) =0.0 
H£ X/ (1+6(X) )} = 0.0 
P(0/Y) =0.0 


The boundary conditions at exit are the same as those for the 
symmetric channel. 


In this case the X,Y coordinates are transformed to T) / g 
coordinates by the relations 


-77 = X/ g=Y-5(x) (2.4.2) 

where 6(x) = 6 ^(x)/L instead of 5(x) = 2 as earlier* 

Following the same mathematical steps as in Sec. 2*3/ the 
following expressions are obtained with reference to Fig* 3 * 


i) Unit normal vector perpendicular to | -constant lines 


-/3 ^ + ©Y 

“--yvr- 


( 2.4*3) 


ii) 


where 3 = ^ 
and a = 1 + 

Unit normal vector tangent to 


( 2 *4 .4) 

( 2* 4* 5 ) 

the 'i? 7 — constant and g —constant 


lines respectively 





(a) Shape of control volume (b) Shape of control 
in physical space volume element 

(dz = unit ) 


Fig. 3 Shape of control volume for non- 
symmetric channel. 
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= e- 


Y ^ ^ 


7] 


,1/2 


a 


(2*4.6) 


iii)' The velocity components^ in the 77- and | -directions, 
respectively, are 


V. 


_ a 


1/2 


u , Uj = -eu + V (2.4.7) 

iv) The moment™ equations In the n- and J -dlreotions, 
respectively, are 


+ V ^ 


dU: 


du, 


^ dX ^ 3 y " 


S ^ + J. (!% . ^^S) 

ax 5T + Re 


( 2.4.8) 


- u(u i (u ^ + 2 

dx Re 3x2 ^ 


Control volume form of the governing equations 
§ -momentum equation : 


(2.4.9) 


/ (U.n) Ug dS - / (n.vUg) dS =r / (3 II 

- ^ U(u dv t i (u g 4 2 H) av 

7?-momentum equation : 

/ (tJ.n) U dS - i- / (n.vu) dS = - ; |g d¥ 

S s V ^ 

energy equation : 



d¥ 

( 2.4.10) 


( 2.4.11) 


/ (U.n)H dS - X / (n^vH) dS = o 

s _ s 

continuity equation : 

J (n-if) dS = 0 
S 


( 2-4*12) 


( 2.4.13) 



32 


vi) The quantities needed to evaluate the surface and volume 
integral in Eq« (2*4#10) are 

/ 3 9 9 

® ®x + 5f 

volume elfanent d¥ = dr)d| 

For surface 

Unit normal vector n = e^ 
u,n = u , 

^ aufc 3Ug 

= ~~ - 3 and dS = dg 

For surface 

n = (-3 therefore U.n = 

n.VU = - 3 and dS = dV 

s 3 s oV 

For surface and dS is identical to those for and 

with the exception that the outward normal n has the opposite sign* 

With the aid of above expressions^ the first and second 
surface integrals in Eq# (2'*4*10) become 

I CU-n)Ug dS = J UUg dg + / Ug U^ dTh- I UUg dg - I UgUg df? , 

^ ^2 ^3 ^4 

/^)dS=/ (^-3^) d§./ (a^-3^) di, 

S ^1 2 

9Ut 3Ufc 3Ufc 3Ufc 

- 1 ^ as ’ i ‘“as"® 
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and the right side of the same equation becomes 


" V d?7dg - / ^ d77d| + ^ / U ^ an d§ 


31' ^ “ ^^2 -'-s • Re'^ - 


+ -1. J- ^ ^ ^ . 2 r o 3U d^6 

3 77 dn^ ^ “ Re (i ^ 3 § ^^2 


2 

drT 


(2.4.14) 


Due to the same reason as before^ the common form of 
Eqs. (2.4.10), (2.4*11), and (2*4.12) becomes 

S {in <P + Cir-i-^)} dg - ; { (u<jp +P(r+^)} d§ + / {u m 


+ P(q+A) 1 d7?-/ {U.g+p(Q+A)} d77=b 

^4 ^ 

where (p stands for U or or H, 

n = 1/Re for momentum equations 
= 1/Pe for the energy equation, 

and 


(2.4.15) 


= - a ^ 

31 ' 


-,9^ 3^ 


(2.4.16) 


For g-momentum equation b is given by Bq# (2.4.14). For T>- 
momentum equation b is given by 


- / II dV = - / (i| - 3 ||) d-71 dg 




3^ 


( 2 . 4 . 17 ) 


For energy equation b is zero. Finally, the continuity equation 
Eq. ( 2 . 4 * 13 ) may be ejcpressed as 


S U d| + / U. dn - / U dg - ; U. d^ = 0 

Q ^ O S Q S ^ 


(2.4.18) 
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The control— volume formulation for the flow through 
channels and pipe is now complete* In the next chapter attention 
will be focussed on the finite-difference approximation of the 
control vol\ime form of the governing equations derived here# 



Chapter 3 


FINITE-DIFFERENCE SOLUTION 
3.*1 Domain -Discreblzat ion : 

For the derivation of finite-difference equation^,, the 
computati-n ai domain in the R - | coordinate system (which is 
rectangular in the 71 - | plane and non-rect angular in the X-Y 
plane) is first discretized^* Only gradients of pressure appear 
in the momentum equations and pressure does not appear 
explicitly in the continuity equation, although it is the 
continuity constraint that is used to determine the pressure* 
Thus if the velocity components and pressure are computed at 
the same location and interpolated linearly for evaluation at 
the control volume faces, the resulting discretization equations 
can permit a physically unrealistic solution [7. ]• This 
possibility of an unrealistic solution can be prevented by 
using different grids for the velocity components and pressure 
and suitably staggering than relative to each others The 
locations where the pressure and other dependent variables 
( except the velocity components) are calculated are designated 
as main grid points^ In practice, for computational domain 
discretization, the main control volume boundaries are first 
drawn in an arbitrarily nonunifoim manner and the main grid 
points are then placed at their geometric centres* In the 
above practice, the main control volume faces are not located 
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midway between adjacent grids points^ which represents less 
accurate finite difference representation of the derivatives# 

In spite of this, the above practice is adopted here due to the 
following desirable features : 

a) the grid point at the geometric centre of the control 
volume represents the control volume better than any other 
point, 

b) half control volume does not occur adjacent to a boundary, 

c) discontinuities in thermophysical properties, boundary 
conditions and the source terms are easier to handle# 

Due to mathomatical singularity, the grid spacing near 
the entrance to the channels or pipe should be small^ Also 
smaller grid spacing is needed near the walls because of larger 
gradient of the flow parameters in the direction normal to the 
wall « 

The grid points corresponding to the velocity components 

are displaced with respect to the main grid points in such a 

way that they are located on the faces of the main control volume* 

The U-control volumes are staggered horizontally (i.#ej., in the 

■T] -direction) with respect to the main control volume, while 

the Ufc-control volumes are staggered vertically (i.#e* in the 

s 

§ -direction) • This staggering is done in such a way that the 
displaced faces pass through the main grid points, while non- 
displaced faces lie along the main control volume faces# Hence 
the resultant shape of the control volume used to ccampute 
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velocity components and other dependent variables becomes L- 
shaped as shown in Fig, 4, In this figure the main grid points 
are denoted by dots identified by letters, P,W, E, S,N/SW, SE,NW, and 
NE whereas U~velocity locations are denoted by arrows (X-») 
identified by P", E", S", SW^', NW" and NE" and -velocity 

locations are denoted by arrows (il) identified by P'/W', E', 
SW',SE^,NW^, and NE'* The U and locations are referred to 

as staggered points.* In the same figure (Fig. 4 ) the main control 
volume, U— control volume, and U. —control volume are, respectively, 
identified by inclined, vertical, and horizontal hatch lines* 

The type of control volume used for velocities near the boundaries 
is different as shown in Pig* 5* A complete picture of the 
discretized computational domain used is shown in Fig* 6 a, in which 
NX and NY are the total no* of main grid points used for 


discretizing the computational donain, respectively, in the 
and § -directions* In this figure boundaries of the main 
control volume are shown as dashed lines* Discretized fojmi 


of physical domain for both symmetric channel or pipe and for 
non— symmetric channel are also shown in Pigs*7a and 7bj, The 
various geometrical quantities needed for discretizing the 
governing equations are also shown in Figs* 4, 5,and 6 in 
FORTRAN variablesj* Attention will now be focussed on deriving 
the discretized form of Eqs* (2*3*55) and (2*4*15)* 

3*2 



The control volume form of the U-moraentum equation is 




SYV(2) , YV(IYPI) 





Fig. 5 Type of control volume at the boundary. 
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Fig. 6(b) Various grid quantities. 



ig. 7 a 


DISCRETIZED PHYSICAI DOMAIN FOR 
A SIMMS TRIG CHANNEL OR PIPE C FOR 
1st CTOIE ONLY) WITHM. = 0.10^ 
p®=1.02, F2=0.98, NX=58, and NY=20 



Pig. 7 b DISCRETIZED PHYSICAL DOMAIN FOR 
A NON-STMMETRIC CHANNEL (FOR 1st 
CYCLE ONDY) WITH 10, P.=1.02, 
P, = 1 .01 , AH3 NY =20 ^ = 
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given by Eq. (2*3.*55) in which cp = U, p= 1/Re and b is given by 
Eq,, (2-»3»56). The control volume used in this derivation is the 
U -control volume such as the one that surrounds the grid point 
P" in Fig* 8* The eight neighbouring grid points of P" are 
denoted by W'', S'', S'*, N", SW", SE", NW" , and NE"* The basic 
principle in the discretization of the surface integrals in 
Eq.* (2«3»55) is to approximate U and its normal derivatives at 
points E^ n", P, and s" on each face between adjacent u— velocity 
grid points# Then the integration is performed by regarding those 
values as constant over the respective faces* With this 
approximation the surface integrals in Eq# (2»3*55) reduce to 

[CUU6 + r (aS)P - (WB + ^ r)p 

A5u + + r a) - (UgU + Is a) ,, 

II m s 

+ - <A (3-2.1) 

where = / ||| (6g)^"^6 dg d7? 

^ U u s 

- 6 d| dn (3..2*2) 

and Q=-ff|/r = ~6ff/A = pP,f=|3|| ( 3-2-3) 

The choice of particular scheme for discretization of the left 
side of Eq* C3#2#l) depends on the relative importance of tbe 



sw Tsw'' s js" 
i I I 

Fig. 10 Main control volume 
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convection and diffusion terras# When convection is small^ the 
central difference scheme can be used to yield a result of high 
accuracy for a suitably small mesh size# However, for fluid 
flow problans/ in general, convection is usually large,# Therefore, 
the scheme should account for the special influence of the 
upstream values# In the presence of convection the particular 
discretization technique used here draws heavily from the 
material set forth in [ 7 ] • 

It is convenient to define the total (convection plus 
diffusion) flux J on the surfaces and as 

Je = (UU6 + i rig (6S)E'^i5^ 

6 3Un _ 

= (WJ6 

= F-g Ug — Djg(Ug,, *“ Up„) (3#2#4) 

and 

= ‘V sfV 

= ^n^ "" ^n'* *” (3*2^5) 

where Fg and P^^,, are the convective fluxes and Dg and D^„ are 
the diffusive fluxes and are given by 

Convective flux : 

Pg = (U6 )e (5§)g’'^Ag^ = (U6 )e SYG(IY) 


( 3-2-6a) 
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SXUClx) 


( 3^2*6b) 


Diffusive flux 




) 


K -1 < 5 e SYG(IY) 
® DJOJdXPl) 


D 


(r)^// (6£) 


n' 


Re '6 n 


K -1 

n" 


sxltC IX) 
DYG( lYPl) 


(3.2-7) 


where SYG(IY)^ SXUdx), DX[JdXPl),and DYGCiyPI) are the FORTRAN 
variables and are explained in Figs. 4. Similarly other FORTRAN 


variables will be introduced later. 


In a similar way the total flux on faces and is 
defined as 


Jp = Fp Up - Dp (Up^ - 

and ~ ^s"^' *” 

where convective flux : 

Fp = (U 6 )p (aDp-^ SYGdY) 

( 3- 2.8) 

^s" = sxudx) 

and diffusive flux : 

6 p SYGdY) 

~ Re ^^^^P DXUdx) . 

( 3- 2-9 ) 

-1 .as . ..sK-1 sxudx) 

^s^ ^6 dygCiy) 
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With this definition, Eq* (3«2*1) becomes 

~ “■ ~ ^ (3.1* 2 *10) 

where S is defined by 

S =[(i*)p (66 )K-i - (X*)g (as)K-i] 

* (6 1 )^ 7 } - (X A)„„ (6 S)^ZW^’lu 

= If’p ■ <^lrh syG(iY) 

+ SXU(IX) 

( 3*2.11) 

Also, from the continuity Eq. (2.3.57), we have 

(U6)e SYG(IY) + (Ug)^„ (6 5)^^^ S}{U(IX) 

- (U5)p (6§)p’"^ SYG(IY) - SXU(IX) = 0 

5 S S 

Using Eqs. (3. 2..6) and (3. 2..8), the above equation becomes 

^E + ^n'' “ *" ^s" = ° (3.2.12) 

With this, Eq* (3.2.10) can be rewritten as 

( Jj, - Eg Up,,) + (J^„ - Up,,) - (Jp - Pp Up,,) 

“* "" ^3** =5 S + ( 3.2. 1 3 ) 


It is shown in [ 7 ] that for any J ( say Jg) 
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( 3-2*14) 


where a^,, = % A j ^ 


El 


1 + IE- Fe,o j 


where Aj^j depends upon the typ^ of scheme considered^ e^g#/ 


for power law schane A | 
for upwind scheme - A | 


D 


El 

'e' 


= 07° ' (1 - 0-1 l^i>^ J 

E 


D 


Ei 

e' 


( 3 • 2* 15) 


= 1 


The notation Ea,bJ means that the larger of a and b is to be 
used* The final form of discretized U-raomentum equation may 
be written by substituting for Upw) from Eg* (3-2-*14) 

into Eq* (3-2-13)- This yields 


- a 


or ap,^ Up„ 


where a 






aj^-r = 


-^E' 
;> , 

^pl 

,F 


n' 







CUp„ 

“ Uj^//) ” a^^(U^ — 

Up^) 



(U3„ 

- Up^) = S + b^ 





+ S + b 
u 





(3-2. 

»16 ) 

IT 

- Pg 

. ol 



E 

^p / 

on 



-f 

E- 

^ 0 n 

( 3- 2< 

.16 a) 

+ 


. 0 J 
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Again can be split up in the following form assuming the 
value at the grid point to prevail over the control volxime 
surrounding it 

= (61)^:^ (Pp - Pe) «p„ f If ( 6 f 6 as an 

Using the above form of b^ ^ Eq« (3»2»16) reduces to 

9-p// Upyj, = 3-g#/ t 

+ {( 6 g)p«^ 6 p« SYG(IY)} (Pp - Pg) + Sb (3,2-17) 

where Sb = S + / ~ ( 6 |)^ ^ 5 d§ d^ (3-2-18) 

y; 5 o% 

and where S is given by Eg- (3-2-11)- 

The last two terms in the above B 5 - (3-2-17) are known as source 
terms- The term Sb appears due to the non -orthogonal nature 
of the T) - I coordinate system - 

For discretization of the source term Sb linear 
interpolation and finite-difference formulae for the variable 
and its derivatives which are not directly available at the 
location being considered are used- 

Using this^ the discretization form of Sb becomes 

Sb = SUW - SUE + SUS - SUN + SUP (3-2-19) 

where 

SUW = (i B ||)p (sg)f-^ SJSUY) 

1 j eSH') 4. (^) \ ( 6 1 )^"^ syg(iy) 
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_ 1 
" 2 


SUE 


SUS 


SUN 


3p 

Re 





U 


SW" 


] 


(5g)p“^ SYG(IY) 
Y(IYPl) - Y(lYMl) 


( 3^2.20) 


‘li If’e ‘5S>e‘^ SIG(IY) 


1 


= 1 Ri c ‘lr>p" 


.3U 


xK-l 


1 


2 Re 


[\-,u - Uq/, + U 


(6§)^”^ SYG(IY) 
^ JJ 1 

e ■}?/# J 

Y( lYPl) - Y(IYNIl) 


NE" "^SE" 


( 3*2,* 21) 




r /^U^ . , /3U 


r 4- f f^) - C^’) 1 

Re L S" P" ^dV S" ^ 

(6t)^7/ SYU(IX) 


( 1/2)SYG(IYM1) 
DYG(IY) 


1 


0 gjf.» 

R^ L^SE^ 


Tuc 


^SW" ■*■ ^^E" “ ^W" “ ^SE" 


(6t)^Z^ SYU(IX) 

O 

XQ(IXPI) - XU(IXMI) 


(1/2)sYG(IYM 1) 
DYG( lY) 

( 3*2*22) 


] 


|w*n“ 


n' 


3 


_nfr(^) .C(^) 

Re p" dyg(iypi) 


SXJ(IX) 


n 


-fe-' 


u. 


Re ^W" 


^^NE*' “ ^NW" ■" ^E" 


(V^)™ ^ 


sxu( ix) 


( 3*2^23) 


DYG( lYPl) 


XU(IXPI) - XU(IYMI) 
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and 


SUP = ; £ (||) 6 di dv 

V 6 y ^ 


<5p» 3l 


^7/ 5n„ SYG(IY) SYU(IX) 


^ , ..^P^ .^P^ ,(V2)sxg(ix)^ 

6p// a§ ag 35 p sxu(ix) 

6^,, SYG(IY) SXU(IX) 


. r (1/2)SXG(IX) 

— [ P„ - Pg + (PjjE - Pgg - P^ + Pg) ] 


K -1 


( 6 < 5 -n/» 


S ^ ""NE ” ^SE "N 

syg(iy) sxu(ix) 

Y(IYPI) - Y(IYMI) 


( 3-0 2^ 24 ) 


The total source term of U-mom^tum equation is 


= Sb + {( 6 g)p“^ 5p,, SYG(IY)} (Pp - Pg) (3-.2^25) 

where Sb is given by Eq,. ( 3*2^19 )• 

The source term can be expressed in the general 
form [.7 ] 



^u “ ^c"'^ ^P" ^P*^ 



( 3»2#26) 

where 

S^„ = Sb + { (6 £ )pZ^ 

6 p#^ SYG(IY)} 

(Pp - Pp) 

(3#2^27) 

and 

Sp## = 0,0 



(3# 2*28) 

With these# Bq« (3£*2^17) can 

be written as 
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= D , 

e' 

^ IdT 

■I + f -''e' 

^W' 

= D / 

w' 

w 

1 ^ tv ' 


- % * 1^1 ■ 

IT “% . 0 

^S' 

= DpA 

l^p 1 

Ie;I - 

T^p , 0]] 

^p/ 

= 

+ + a^y + ag_, 


where ^e' ' V / 3.nc3. 

are the convective and diffusive 

Convective fluxes t 


/ 0 T} 

ol 

J] (3.2*32) 

' respectively 

fluxes and are given by 


^e' 

=. (U6)g/ (a5>e7^ 

SYV( lY) 

^w' 

= (ua)^ (as)^"^ 

SYV( lY) 


= (Uj)n (ag)^"i 

SXG( IX) 

% 

= (Uj)p (a5)f"^ 

SXG(IX) 


where the velocity j 

the following relation 


( 3- 2^33) 


and (U^ )p 


are obtained from 


U . = + (Uvr,^ - 


(1/2)SYG(IY) 


p.- 'pjw T '‘‘-'N/# ^pjw^ , . 

e ^ ^ DYG(lYPl) 


(1/2)SYG(IY) 

W wwr DYG(lYPl) 


(Uj)^ =(V2)£(Uj)p, + 
(U|)p =(V2)f('^j)s' + 


( 3.2.34) 


( 3.2f.35) 
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and diffusive fluxes 


V = ^ 


,K-1 


,, SYVdY 
DXG( IXPl) 


w' Re v/ 

1 ^ R 1 SXG(IX) 

% = :& 

<5 ^ DYV(IYP: 

SXG( IX) 

'S = ie ^ ^ 

i<e 6 ^ DYV(iy) 


D 


Eq* (3»2*30) can also be expressed as 


( 3-2*36) 


^U| - 3 -e' ■*" '*■ 

+ asdUg)s^ + (6g)pr^ 5p, SX(^(IX)} (Pp-P^)+Sb 

(3-2^37) 

where 

Sb = s + ; 3 |l (6§)^''^6 d§ dh - / [u(^ il + 

V ^ ' \L O O.IJ 

( 6 ae S’! + ^ [u(S 0 + ^ (g^)^ ' ^ 0 

+ at’3 (6S)''"^6 d{ 377 + ^ ^ (g ^ - I |2) 

^If 

- h i ^ "i e“>] i6S)^-^6 dS dv 


( 3# 2.38) 
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The quantity S in equation (3.2.,:^) is given by Eq. (3-2.31)- 

Using the same concept that was used for discretizing 
the source term in the U-momentum equation/ the discretization 
form of Sb is given by 


Sb = SVW -- SVE+ SVS - SVN + SVPl 

- SPRIME^ - SVP2 + SVP3 +SVP4 

+ SVP5 - SPRIME^ (Ug)p/ - SVP6 (3.2-39) 

where each term on the left side of Eq, (3-2-39) is as follow: 


au 


^ SYVdY) 


3, . aUfc au. au. 

— 2L. r ( — s.) , + f( — M,^,) 


~ SXG(i:XMl) ^ ^ . 

^ •] (<5|)J;7^ SYV(IY) 


DXG(IX) 


J SXG(lXMl) , SYV(IY) 


+ (Ut)oT.T/} 


K-1 


^ DX3( IX) 




^ YV(IYPI) - YV(imi) 


SVE 


au 


= (fea^)e^ (6^)^7^ SYVdY) 


au. 


au. 


SXG( IX ) 




Re '-'as 
x(6S)^7^ syv(iy) 



+ 


svs 


SVN 


SVPl 


■i SXG(IX) ^ . SYV(IY) 

(U.)g,} ^ — (3.2^41) 

DXG(lXPl) ® YV(IYPI) - YV(IYMI) 

O 9Ufc ^ ^ 

= (&3Tr>p 

3p t r/- M 

= Bil f<37r>P' + SX3(IX) 

” "5 Re ■*■ ^SE' " ^SW' ] 

. SXG(IX) 

(6§)p (3#2#42) 

X(IXPI) - X(IXMI) 

3U. 

= <teFir>N 

= 1 S 

" "i ^W' **■ ^NE' “ 

SXG(IX) 

X(6g);Y'‘^ (3.2*43) 

X( IXPl) - XClXMl) 

= ^ 3 1^ (6^)^*"^ 6 d§ dT? 


= 3p/ (6S)p7^ 6p, SXG(IX) SW(IY) 



f ( 


~ (ls)p} 


dv' 


1 syg(iy) 

DYG(IYPI) 


6p^ SXG(IX) SYV(IY) 


X(6S)^7^ 
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6p, SXG(IX) SYV(iy) 

X -Z 

X(lXPl) - X(IXMI) 


‘NW 




•| sygCiy) 
DYGClY^’i 



K -1 

P' 


( 3*2^44) 


SPRnME^.(Ug)p^ 

t^P* 

= Up^ -r^ 

u p# 

. SPRIME- = U„. 

» # J. ir 


■ ^m'^P' ( 5^)p7^ 6p/ SXG(IX) SYV(IY) 

“ (|~)p/ (6g)pr^ 6p, SXG(IX) SYV(IY) 

(3-2*45) 


2 

SVP2 = S U(Ug (6g)^~^ 6 dg dT] 

V . dT7 


2 

= U^, gp, ^~“|^P' (ag)^"^ SXG(IX) SYV(IY) (3-2*46) 

d 7 ) • 


'p/ Up, 

The Up, in Eqs- (3- 2-44) and (3-2-45) is given by 


U 


P' 


1 SYG(IY) 

= [Up + (Un - Up) ^ ^ ] 

DYG(lYPl) 


"I [Up« + u^,, + (Ujj,, + Uj^^ 


^ syg(iy) 

- U„„-U„^) 2 ^ -__] 


P" w* 


DYG(lYPl) 


( 3-2.47) 


= :& 4 ['J + fl - f i 

4 + (K— l) — (6$)^ ^5 dg dTj 

dri^ ^6^ 


u 


P' 


Re 


{^p. 


(di) 

an 3 



d 6 ^ 

3 t)^P 


3 gp/ 

6 p# 
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(K“i) - — ^ o c^)p,} 


dr/ 


2'P 


6p/ 


dn P'' 


^(6g)pr^ 6p/ SXG(IX) SYV(IY) 


( 3.2*48) 


where Up, is given by Eq, (3*2.47) 


f 6 d§ d77 


SVP4 = ~ / r ^ 

- _2_ (9U) r. ^d^a) _ rd6) 1 

- Re "V (Sl)p^ 


5p, SXG(IX) 


where (|^)p# 


dv^ 

SYV( lY) 

f(lH) . (SLE) , I 
an N an dtc(iypi) 


(§ii) + f (IH 

Sr?^P ^ 




i syg(iy) 

j 

dyg(iypi) 


SXG( IX) 


( 3.2.49) 


(3. 2.. 50) 


SVP5 = 


Re 




= ^(IS) 


a 






_£. (~) f—i— f^) - ^ P'^ ( ^^^ ) 1 




X ( 6§) pr^ <5 D/ SXG( ix) SYV( IY) 


(3.2.51) 


where (g^)p 4 = ^ ” ^P" ^NW* 


Uw^) . 


" DyG( lYPl ) 

K-1 


SPRIME2(Ug)p^ = (j ^ ^ (6S)^”^6 d§ dT? 


( 3 # 2#5 2 ) 


= i 'K-D p-i^ (Uj)p, (6S)fr^ 6j, 

Sp# Op/ 


SXG(IX) SYV(IY) ( 3 . 2 . 53 ) 
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,•*. SPRIME- = (K-l) 4- 
^ R© 


(Sp. 6p,) 


2 (6^)p7^ 6p/ SXG(IX) SYVdY) 

( 3.2-.54) 


SVP6 


= ^5^- {(K-l) ^4^ I3U} (6^)^"^ 6 dt 


(K-l) 


:252 


Re 


(^ p/ '^p' ^ 


2 Up, G6§)pd 6p^ 


K-l ^ SXG(IX) SW(IY) 


/# The total source term is 


= Sb + {~ (6S)p7^ 6p, SXG(IX)} (Pj 




In generalized formas can be expressed as 

U fc 


= ®C' V 


( 3.2*55) 


( 3 #2*56 ) 


(3.2.57) 


where S^, = SVW - SVE + SVS - SVN + SVPl - SVP2 


+ SVP3 + SVP4 + SVP5 - SVP6 + { 


SXG(IX)} (Pp - Pj^) 


f cfcNK-1 

677 ^P' 


and Sp, = -(SPRIME^ + SPRME^) 

With these, Eqj, (3i*2-37) may be written as 


(3.2.58) 

( 3.2.59 ) 


( arj, — S-p, ) (Ufc)-r)/ = arj,, (Uf.)-R>, + 


(U*,), 


a.p# ^ V'-'g^p/ = a.£^ T ^ 

+ ( 3.2.60) 

c) Discretization of energy ejuation 

The derivation of the discretized equation for enthalpy 
is based on the assumption that the flow is available a priori. 
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In this case Eq, (3.»2^55) is once again considered, but now <P 
is taken to be H and b = 0» The control volume of interest 
in this case surrounds the main grid point P as d^icted in Pig,# 1C 
The eight neighbouring points are identified by W,E, S,N, SW,SE, 

NW and NE. 

Now using the same mathematical steps that were used for the 
derivation of discretized U-momentum equation, the discretized 
form of energy equation becomes 

= ^'W ^E ^ 3»2#6l) 

where S is given by 


S 


= h&|f>w (8 5)^-1]S«5(IY) 


( 3-2#62) 


and the coefficients a^ , a^, , ag , a^^^and ap 

= “w ^ . 0 1 

w 

^E=“e*l^! +ir-I'e - oJ 

G 

= I>s •^1^1 + ' °3 

s 

^Id^I + E-^n ' °3 

n 

and ap = a^ + a^ + aj^ + Sg . 


are given by 


( 3-» 2«6 3) 


The convective fluxes P^ / ^ ^nd P^ and the diffusive 

fluxes Hi given by 
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Convextive flux : 


^e = 

(W>e 


SYG( IY) 

= Up,, 6^ 

e 

= 


(64)^-l 

SYG(IY) 


(6£)^-' 

= 


( 

SXG( IX) 

= (Ug )3, 

(6£)^- 

= 


(6£)K-i 

SXG(IX) 


(61)^- 


^ syG(lY) 
SYG( IY) 


Diffusive flux : 


1 TC-I SYG(IY) 

£) — 1 ( 6 £ ) ^ w 

w Pe w dX3(IX) 

6 SYG(IY) 
dxG(IXPI) 

DYG( IY) 


D. =4 (?)„ (6iL‘^-i 


"n = ^ 6 'n 


n 


DYG(lYPl) 

The discretized form of source term S is giver by 
S = SHW - SHE + SHS - SHN 

where 


SHW = ||)^^ sro(iY) 

= If Kw-«sw+ '%-Hs'%W + Hsw) 


( 3*2-64) 


( 3— 2-#65) 


( 3* 2-66) 


•| SXG(I}M1) 
. DXG( IX) 


] 


SYG( IY) 


X ( 6 S)^"^ 


Y(IYPI) « Y(IIMI) 


(3.2*67) 
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SHE = (^1^)^ SYG(IY) 




^ SXG(IX) 

H - + H ) 2 : ] 

^ DXG(lXPl) 


. SYG(IY) 

X ( 

® Y(lYPl) - yCIXMI) 


( 3.2.68) 


SHS = (^- 15) SXG(IX) 

Pe dV s ^ s 


= [H 


“ ^TaT "* T + ^^OT.t) 


1 sygCiy^i) 


“ Pe L“SE “SW “SE DYG(IY) 

„ . SXG( IX) 

X (5 § ) _ 

^ X( IXPl) - x(lXMl) 


SHN = (r§- ^) (6?)^’"^ SXG(IX) 

Pe a 77 


p| [He - % + (h^e - Hhw - He + h,,) 

SXG( IX ) 

X(lXPl) - X(lXll) 


The total source term in this case is 


Sr = S 


Again the general form of Sj^ is given by 


■| SYG(IY) 
DYG(lYPl) 


] C6§) 


(3.2.70) 


wh.ere 


% = ®c ■*■ ^P ^P 


S and So = 0 


( 3.2.71) 
( 3-. 2 .7 2 ) 


With these, Eq. (3.2.61) may be e3q>ressed as 


(ap - Sp) Hp = ag Hg + ^S ^S **" ^ 


( 3 . 2 . 73 ) 
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d) Discretiza-blon of continuity equation i 

The control volume of specific interest in this case 
is that which surrounds the main grid point P in Fig. lo* The 
surface integrals in Bq« (2»3»57) are now approximated by- 
evaluating the integrand at points e ^ n ^ w , and s respectively* 
Then the integration is performed by regarding those values 
as constant over each face to obtain 

SYG(IY) + (Ufc)„ SXG(IX) 

- u Ar SYG(IY) - (u-)^ SXG(IX) = 0 

WWW 500 

or 

Up« SYG(IY) + (U^)p, ^5^^n ^ SXG(IX) 

- 6^ SYG(IY) - (Ug)s. (<5g)^r^ SXG(IX) = 0 

(3.2#74) 

Using Eq. { 3f2fS4) ^ the above Eq. (3»2*74) may be written as 

Pe + Pn - ^w - Eg = O (3^2.^74a) 

Before deriving the pressure-correction equation at this stage 
let us discuss the necessity of using under-relaxation factor 
for solving Eqs* ( 3*2*29)/ (3*2*60) and (3*2*7 3)* 

e) Under-relaxation : 

The generalized form of Eqs* (3*2*20), (3*2*60) and 
( 3 * 2 * 73 ) for a control volume is given by 


(ai - Sj_) ^ ^nb ^nb + ^<^1 


(3*2*75) 


where stands for U , , or H , 

or Sg and 1 indicates the point P-' 
for the U-inomentum, the U|-inomentum/ 


stands for 

^ P' or P , respectively, 
or the energy equation# 


Equation (3»2*75) appears to be linear# However, the 
coefficients of Eq# (3#2#75) may themselves depend on one or 
more of the dependent variables represented by cp # To account 
for the resulting inter- equation linkages and nonlinearities, 
repeated solutions of the nominally linear form of Eq# (3#2#75) 
are required^ Each of these iterative solutions is defined 
herein as a ^ cycle*'# At the beginning of each cycle, the 
coefficients are evaluated using the cp values obtained in the 
previous cycle# With the cycle— by-cycle change in coefficients 
of Eq, (3#2^75), the resulting changes in the (P values can be 
quite large, and this may cause slow convergence or even 
divergence.# To moderate the changes in successive solutions 
for (p, and thereby improve convergence, under— relaxation is used# 


Patankar [ "7 ] introduces under-relaxation into Eq# (3#2#75) 
through as follows i 


( a^ — ) 


1-a 


a 


R 


(Pi = 


^ \ib ^nb 


R 




a. 


(a. 


R 


sp 


( 3.2*76) 


where <Pi is the value of ^i from the previous cycle, whereas 
Van Doormaal and Raithby [l2] introduce under-relaxation through 
use of an E-factor according to the following revised form of 
Eq # (3»2#75)# 
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(ai-Si)(l + |) <P^ = S + s + I cp; 




or D^ <p^ = S 0 )^^ + ®<|.^ + i (»1-®1> fl 

where D^ = 

From Eqs* (3*2*76) and (3*2*77) we may write 


(3.2.77) 

(3.2.78) 


Under— relaxation through the use of E-factor has better physical 
meaning than that using <Xj^ ^ as discussed by Van Doormaal and 
Raithby [12]. In order to accelerate convergence/ values of 
E well in excess of unity are desirable. In fact values of 
E in the range 1 to 20 were found to be useful for the present 
problems. 

f ) Derivation of pressure-correction equation 

We use the SIMPLEC procedure [l2j for handling the 
velocity-pressure linkages. In this method the pressure field 
is first guessed. With this guessed pressure field/ coefficients 
of the momentxam equations can be evaluated allowing these 
equations to be solved to obtain the flow fieldj* In general 
this flow field does not satisfy the continuity equation ( 3*24*74) ♦ 
Therefore/ this guessed pressure field is corrected so that 
the resulting velocity field satisfies the continuity equation. 
This is accomplished by the pressure correction equation which 
is derived in accordance with [l2] by combining the continuity 
equation with truncated forms of the momentum equation. After 
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solving the pressure correction equation following the 
recommendation given in [12] y the velocity and pressure fields 
are correctedy and the procedure is repeated until the flow 
field satisfies both the continuity and momentim equations* 
Details for the above outline fellow. 


Let us say that for the guessed pressure distribution P* 
the U'^ and velocity distribution obtained by solving the 
U-momentxM and Ug-momentum equations respectively satisfies 


Dpy* = s ^ SYG(iy)}(p5-P*) 


Vl 

i (a 


P" 


P^/ 


,) u; 


‘-'p/y 


■P E' 

(3^2 # 800 . ) 


and 


oc 

Dp,(Upp, = 

( Pp - ^ ( -ap/ “ 

where Sk^ and Sb^ are given by Eqs* 
respectively. 


(6g)p7^ 6p/ SXG(IX)} 

Sp/)(Ug)p, (3-2-80b) 

(3.2.19) and ( 3 .. 2* 39) 


while the U- and -velocities obtained from Eq. (3.2.75) 
using the correct (but unlcnown) pressure distribution P satisfy 
the continuity equationy the U — and U^— velocities obtained 
from Eqs. (3*2*80) do not in general satisfy the continuity 
equation. Correction of the guessed pressure by P' = P - 
is therefore necessary to correct the U field by U'' = U - U 
and field by U| = Ug - U*. The relation between P* and U' 
is obtained by subtraction of Eq, (3.2#80a) from Eq. (3*2.75) 


leading to 
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Dp„ = E + {(61)^,';^ 6p,, SYG(IY)}(P'-Pp (3.2-81) 

The pressure P and velocity U that satisfy both the continuity 
and the momentum equations are 



( 3- 2-8 2a) 

p. 

( 3-2-8 2b) 


Attention will now be given to the method used to find P'^. 

The exact equation for P' derived from Eqs. (3-2-81) 
and (3- 2-8 2a) and the continuity constraint, is complicated and 
unsuitable for economic calculation- In the SIMPLE method [7] 
the term E ^nb (3-2-81) is ignored whereas a term 

of similar magnitude on the left side of the equation is 
retained, leading to inconsistency of the method^ 

For a "consistent" approximation, leading to a 
suitable simple expression for P', the term Sa^^^Up,, is 
subtracted from both sides of Eq, (3-2-81)- This yields 

(Dpr- “ Up- = 2 a^^ (U^^ - U'„) 

+ {(6S)p”^ 6p„ SYG(IY)} (PJ - P') (3-2-83) 

In this method termed "SIMPLEC" by Van Dooxmaal and Raithby [l2] 
the term S a^^ neglected- With this approximation 

Eq« (3-2-8 3) becones 

(Dp^ ^P" = 6py SYG(IY)} (E^ - P^) 
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or 


11 

pj 

1 

Pi) 

E 



where 



<5p« 

SYG( lY)} 

( 3*2*84) 


Pp,/ - s 

^nb 


Thus, 

Up„ 


(P^ 

- P-' ) 

E^ 

( 3*2.85) 


Eq* (3*2#a5) is an intermediate st^ to link velocity with 
pressure, and is called the velocity correction equation# A 
similar expression can be obtained for 

-- Pp (3.2^86) 

where d^^ = { (6 6^,, syG(lY)}/(D^„ - S (3.2#86a) 

Now in a manner similar to that used for deriving Eqs# (3*2#85) 
and (3*2*86), we have (Using Eq* (3#2.#80b)) 


(Uj 

)p, — (U"g)pr + dp, (Pp “ Pjj) 

( 3-2*87) 

(u 

)g, = "" 

(3-2*88) 

where dp. 

and dg, are given by 



{ (|)p. (6^)p7^ 6p^ SXG(IX)} 

( 3*2*89) 


D-n, — V a , 

P' ^ nb 


{(|)g/ (6§)sr^ 6g, SXG(IX)} 

( 3.2*90) 

o-S^ 

°S- - S ^nb 


Substituting Eqs* (3*2*85) - (3*2.88) into the continuity 


Eq* ( 3 * 2*7 4 ) , we have 



68 


6pj, SYG(iy) + i6^)p7} SYGCIY) dp^ 

- U*,, 6^^ SYG(IY) - 6^., (6 0)^$:^ SYG(IY) ( P^ - P') 

+ (Ug*)p, (6§)p7^ SXG(IX) + (6g)p7^ SXG(IX) dp, ( P' - P^) - 

(Upg, (6£)g7^ SXG(IX) - (6g)s7^ SXG(IX) dg, ( P' - Pp = 0 
from which we can write 

ap P' = ag P^ + P^ + P^ + ag + B (3*2.91) 

where 

= <ap,, 6p„ (5g)p7^ SYG(IY) 

= ’^W' -5^" SYGCIY) 

% = cip/ (6g)p7^ SXG(IX) (3. 2.9 2) 

ag = dg, (6g)g7^ SXG(IX) 

®P = ■*■ 

B = - (u|)p, (6g)p7^ SXG(IX) + (Upg, (6§)g7^ SXG(IX) 

- 6p„ (6S)p7^ SYG(IY) + SYG(IY) ~ 

( 3.2.9 3) 

and where dp„ , d-^„ , dp, and dg, are given by Eqs. (3*2.84), 

( 3.2.86a), (3.2.89), and (3.2.90) respectively. 

The direct solution of Eq, (3.2.91) for all the P* values 
on the line IX in Fig. 6 can be obtained with one application of 
the Thomas algoritlim by suitably guessing the off-line dependent 



variable values^ Such a line-by-line solution is the basis of 
an iteration scheme that solves along each IX-line and then along 
each lY-line, and repeats the pattern until convergence is 
achieved. The rate of convergence of such a scheme depends 
crucially on the treatment of the off line values of the 
dependent variable^ 

Let us suppose that in the current iteration Eg. (3.2.91) 
is to be solved along an IX line^ sweeping in the direction 
of increasing IX, Now according to Patankar' s suggestion [7 ] 
the available estimate of P^ is frcm the previous iteration 
i.e. [P^]° / taut according to Van Dooxmaal and Raithby [l2] the 
best estimate of P^ is 

Pe = ■ LeP) 

* [pj.p + (e - 1) (p' - [p']°) 

where 0 is a relaxation parameter such that for 0 = P^ is 
taken as [p^]®.. The optimal value of 0 depends upon the 
problem.. In general^. 1 < 0 < 2. The best available estimate 
of P^ is that obtained from the just completed solution along 
the (IX - 1) line. It has been observed that Van Doormaal and 
Raithby' s approximation for off-line values accelerates the rate 
of convergence. Therefore, with Van Doomaal and Raithby' s 
approximation Eq. (3*2-»9l) takes the form 

( ap - (0-1) ap) ^ = 9-s ^ + ^N ^N % 

+ ap {[P^]° - (0-1) + B (3.2*94) 

A similar estimation is made for solution along an lY— line. 
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3.* 3 Discretized Equations for t-he Non -Symmetric Channel 

Following the same procedure as in Sec. 3*2/ the 
following equations are obtained for the non-symmetric channels 

a ) Discretization of Utnomentum equation 

The U-momentum equation (2.4.15) with <p = U ] 

can be written in discretized form as 




+ ag„ Ug^ + {SYG(IY)} (Pp-Pg^ + Sb 


ap 


where Sb = S + / /3 dr/ d§ 


( 3.3*1) 
( 3.3.2) 


® If'e] SW3(IY) 


+ fi’s" ■ sxudx) 


( 3.3.3) 


‘Re 37/ n^ 

The coefficients in Eq, ( 3 . 3.1) are given by Eq.* (3.2^16a] 
where the convective and diffusive fluxes are e 5 <pressed as 

Convective flux : 


Pg = Ug SYG( lY) , Fp = Up SYG( IY) 

F^,^ = SXU(IX), SXU(IX) 

Diffusive flux : 

1 SYG(IY) ^ 1 SYG(_IY) 

E = DXJ( IXPI) ' P “ Re DXU( IX)’ 


(■3.. 3.4) 


D 


n 


/# 


^ X a^„ SXU(IX)/DYG(IYP1) 


(3.3.5) 


V = Re °^S" SXU(IX)/DYG(IY) 
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Using the same assumptions as in Sec« 3# 2, the discretized form 
of Sb becomes 


Sb = SUW - SUE + SUS - SUN + SUP 
where 

~ 2P^ *■ "^S* ^NW" " ^SW"^ 


SYGC IY) 


SUE 


1 PE 


" I Re ^’^N" “ ^S'' ^NE^ “■ ^SE"^ 


Y( IYPl)-Y( lYMl) 
SYG(IY) 

Y( IYPa)-Y( lYYl) 


( 3^3.6) 


( 3-3.7) 


( 3.3.8) 


^ [Use., ~ Ug^, + (Ug,, - u^,, - Ugg,, 


■*■ ^SW"^ 


ISYG(IYMI)^ SYU(IX) 


DYG(IY) XU(lXPl) - XU(IXMI) 


(3.3.9) 


“ Re [^E" “ ^W" ■*■ ^^NE" ”■ ^ 

SXU( IX) 


NW" "^E'^ 


■i SYG(IY) 


DYG(lYPl) XU( IXPl)-XU( IXyil) 


( 3.3.10) 


SUP = [%~Ps'^^^ne“^se''^n‘^^s^ sxu( ix) - 


J SXG(IX) syg(iy) sxu(ix) 


Y( IYP1)-Y( lYMl) 

( 3.3.11) 


-*• The complete source term in the U-momentum equation is 
S^ = Sb + {SYG(IY)} (Pp“PE^ 

“ s^ .. + Sp., Up.. ( 3 . 3 . 12 ) 

where S^.. = Sb + SYG( IY) ( Pp-Pg) (3.3.13) 

= U 


( 3.3.14) 
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With these^ Eq, (3»3;*l) may be expressed as 


^ ^p/f — ^W" '^W-" ^N'^' ^S'*' 

(3:»3#15) 

b) Discretization of Ug -momentum equation : 


The Ug-momentum' equation [Eq, (2«4#25) with ^ = Ug j 
can be discretized as 


^E' '*’ "*’ 

+ ag, (Ug)g, + {ttp, SXG(IX)} (Pp-Pjg) + Sb 


( 3:- 3.16) 


v/here Sb is given by 


Sb = s + ; B — drids-zu^ dr) d g 

V 3r) \s ar)^ 


with 




dr? dg+ ^ ^ ^ dr? dg 

dr^ Re^3r?^^2 




Re 


d-r?*^ 


= [{-a.ii). - (^-^> J 

L ^Re a I 'v/ ^Re 9^ 

+ - 'fe ] S>!S( IX) 


( 3.. 3.17) 


( 3-3.18) 


•ReaT^'P 'Rear? 

The coefficients ag, , a.^, , , a^, and ap, in Eg. (3.3.16) 

are given by Eq, (3.2.32) where the convective and diffusive 


fluxes are : 
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Convective flux : 


V = V syv(iy) , SYV(IY) 

( 3.3*19) 

Pn = SXG(IX), Pp = (U^)p SXG(IX) 

Where ^ , (Ug )j^ /and (l^ )p are obtained from Eqs* (3*2*34) 

and (3*2*35). 


Diffusive flux J 


SYV( lY) 

D _ 1 D , = 

® DXG(lXPl) ^ 


^ SYV( lY) 
DXG( IX) 


% = 


“n 


SXG( IX) 
DYV( lYPl ) 


/ — 




SXG( IX) 
DYV( lY) 


( 3*3*20) 


Discretization foxm of Sb : 


Sb = SVW - SVE + SVS - SVN + SVPl - SVP2 

+ SVP3 + SVP4 - SVP5 (3..3*21) 

where 

SVW = [ (Ug)j^^ - (Ug )g^ + { (Ug )j^/ - (Ug )g, 

i SXG(IXMI) SYV(IY) 

^Ug^NW^ + DXG(IX) ^ YV( IYPl)-YV( lYMl) 

( 3*3*22) 

SVE = [(Ug)jj^ - (Ug)s/ + (^Ug)j^g, - (Ug)g£, 

1 SX3( IX) SYV( lY) 

- (Ufc)vT/ + (Ut)c*} ] 

^ ^ S S SXLF(IX) YV(IYPl)-YV(imi) 

(3* 3*23) 
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SVS = - — [(Ug)g, ~ ^SE-' *■ 

SXG( IX ) 

X(IXPI) - X(IXMl) 

1 

“21^ -* + (U^)j;jg, - ^ 

sxgCix) 

X(IXPl) - X(lXMl) 

1 syg(iy) 

SVPI = [Pe-Pw^-(Pne-Pnw-Pe+Pw5 dy^u y ; ! ) ’^ 


(3.3^24) 


( 3.- 3- 25) 


SXG(IX) SYV(IY) 
x(ixPi) - xdxyii) 

SVP2 = u|, (^)p^ SXG(IX) SYV(IY) 
dri^ 


SVP3 



u. 


Pf 


^3. 

SXG(IX) SYV(IY) 

d 77 ^ ^ 


SVP4 


- - 2 _ ( 3 U) 

" Re 


SXG(IX) SYV(IY) 


SVP 5 



(IR) 

P' 


(^) 

d„2 


SXG(IX) SYV(IY) 


( 3*3.26) 

( 3* 3* 27 ) 

( 3*3.28) 
( 3.3.29) 

( 3.3.30) 


where Up, , (|^^p/ / and are given by Eqs. ( 3 .. 2 . 47 ), 

(3.2.50) and (3.2.52) respectively. 


The complete source term in the U^ -momentum equation is 


Su 


§ 


Sb + hp, SXG(IX) (Pp-Pj^) 

S^, + Sp, (u^)p. 


(3. 3.31) 
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where = Sb + cCp* SXG( ix) (Pp-P^^) (3*3-32) 

Sp^ = 0*0 (3*3*33) 

with these, Eq* (3*3*16) may be expressed as 

( 3 p/~Sp, ) (Ug)p, =5 a^i, (Ug )p, + a^^, (XJg)^ 

t a|^j / ( U g ) + ag /(^^^g/ t , (3*3—34:) 

c) Discretization of energy ecfuation : 

The energy equation [Eq* (2*4*15) with = H J can 
be discretized as 

ap^p = ^ (3-3-35) 

where S is given by 


+ - (| 5 ||)„] SXG(IX) 0 . 3.31 

The coefficients a^ , a^ , ag , a^j and ap of Eq» (3-3-35) are 
given by Eq, (3-2-63)., where the convective and diffusive 
fluxes are given by 


Gonv active flux : 

= Up,, SYG(IY) , P^ = U^,, SYG(IY), 

Pg =(Ug)^SXG(lx) , Fjj =(Ug)^SXG(lX) 


( 3— 3— 37 ') 
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Diffusive flux : 

. sygC iy) 


D 


SY3(IY) 


Pe 


DXG( IXPl) 




w Pe 


DXG( IX ) 


D„ = ~ a 
n Pen 


SXG(IX) 

DYG(IYPI) 


,, a 

s Pe s 


SXG( IX ) 

dyg(iy) 


(3.3.38) 


The discretized form of source term S is given by 

S = SHW - SHE + SHS - SHN 
where 

+ '% - Hs - %w + «SW> 

i SXSdXMl) syg(iy) 

X ^ -] 


(3.. 3.39) 


DXG(IX) Y( lYPl ) -Y( lYMl ) 


( 3.3.40) 




4 SXG(IX) SYG(IY) 


] 


DXG ( IXPl ) Y( I YPl ) -Y( I'M 1 ) 

( 3.3.41) 


SHS = ^ [ Hgj, - + (Hj, - H„ - Hgg + Hg„) 

•1 STO(imi) S}il3(lx) 


DYG( iy) 




x( ixpi)--x( ixyii ) 


( 3.3.42) 


0 


SHN = ^ [He - + (HeE - %W - % + 


SYG(IY) 


X 


'] 


SXG( IX) 


DYG ( lYPl ) X( IXPl ) -X( im 1 ) 


( 3.3.43) 


The complete source term in this case is 


Sh = s 


= + SpHp 


( 3.3.44) 
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where = s and Sp = o (3*3»45) 

with these, Eq, (3.3*35) may be written as 


(ap-Sp)Hp = agHp + a^^-^ + agH^ + aj^I^ + (3.3.46) 


d) Discretization of continuity ecruation : 

The discretized form of continuity equation [Eq..(2.4.18) 
is 

SYG(IY) + (Ufc) SXG(IX) - U„ SYG(IY) - (U.) SXG(IX) = 0 
e s Tx w b s 

(3.3.47) 


e) Pressure-correction equation : 

The velocity-correction equations are given by Eqs. 
( 3.2.85), (3.2.86), (3.2.87), and (3.2.88), where 
dg/ ,and dp, are given by 


SYG(IY) 


■^p// 


Dp.r - S a^^ 


SY3(IY) 


( 3.3.48) 


SXG(IX) 


<3-5/ = 


D 




E a 


dp, — cCp, 


, >-jp. 


hb 


SXG( IX) 
Dp , — E a 


nb 


and the pressure correction equation is given by P 
where P'' is obtained from equation 

(ap - (0-1 ) ajg) Pp = agPg + ^ 

+ ag {[Pg]° “ (e-l) [Pp]°} + ® 


= P* + P' 


( 3.3.49) 
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where 

ag = dp,, SYG(IY) , a^ = d^„ SYvj(IY) , 

a^ = dp, SXG(IX) , ag = dg, SXG(IX) (3-3-50) 

ap = ag + a-^ + a^^ + a^ 

and 

B = -(ut)n/ SXG(IX) + (U|)o# SXG(IX) 

- Up„ SY3(IY) + SYG( IY) 


( 3-3-51) 



Chapter 4 


RESULTS AND DISCUSSION 

Results for both developing and fully-developed region 
in a short wavelength sinusoidal symmetric and non-symmetric 
channel and a short wavelength coverging-dlverging pipe are 
presented in terms of computer generated profiles as shown in 
Figs* ll to 71* Figs* 11 to 29 show the behaviour of the 
various quantities for the symmetric channel* Similarly/ 

Figs* 30 to 52 and Figs, 53 to 71 depict the behaviour of the 
same quantities for the pipe and for the non-symmetric channel 
respectively* An earlier studies of flow through such ducts 
are restricted to the fully developed flow in a long wavelength 
channel or pipe of small wall amplitude whereas the present 
analysis is for a. short wavelength channel and pipe with 
sufficiently large amplitude* Also results are obtained for 
both the developing and the fully developed region.* That is 
why the present results are not compared with any pre/ious 
results.* From the nature of the profiles/ however/ it is 
felt that the present numerical procedure can predict the flow 
through such ducts quite accurately* 

The numerical worh was performed for a dimensionless 
wall amplitude LAMBDA (x/l) =0*10, 0*20/ and 0*25* The 
Prandtl number was fixed at 0*72 (for air)/ while the Reynolds 
number/ Re* was given values of 100 and 500# Non-uniform grid 
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spacing was enployed through factors and F 2 as shown in 
the Appendix in both directionse Solutions were obtained for 
a prescribed wall enthalpy only* 

Though it takes about 10 cycles for the flow to get 
fully developed/ the figures depict only the first few cycles 
and the last two cycles* Due to symmetry, as in the case of 
symmetric channel, the figures show the profiles in the upper 
half of the flow domain only* In the figures the solid 
boundary of the channels is drawn using a very small number 
of points* Hence these boundaries have not been r^resented 
exactly in the figures* 

4i*l Symmetric Channel : 

Figs* 11 to 13 show the behaviour of the dimensionless 
velocity components U and V and the velocity vector U for a 
dimensionless wall amplitude X/L = 0*10 and Re = 100,* Similarly 
Figs* 15 to 17/ Figs. 19 to 21, and Figs- 23 to 25 d^ict the 
behaviour of the same quantities for X /L = 0*10 and Re = 500, 
X/L =0*20 and Re = 100, and ,X/L =0*25 and Re = 100 
respectively. It is observed that the negative V-component 
at the minimum cross-section and the positive V-component at 
the maximiom cross-section increases with Re* This is due to 
the effect of inertia [3]* 

It can be easily observed that the velocity ccamponents 
U and V arc affected significantly by the parameter X * It is 
seen from the Figs;* 11 and 13 that for X/D = 0*10 and Re = 100 
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there is no separation of flow* But when Re increases to 400^ 
it was found that the flow begins to separate as ^ at the wall 
tends to zero* Prom Pigs* 15, one finds a small back flow near 
the channel wall for k/L = 0.10 and Re = 500. The co*»»»ei^ ending 
separation point (s) and reattachment point ( R) are also shown 
in Pig* 15* 

A comparison of Figs* 11, 15, 19, and 23 shows that as 
the amplitude X/L increases, the separation point moves 
upstream and the reattachment point moves down-stream. Similar 
effect has also been observed for an increase in Re. Therefore, 
the s^arated flow region grows with wall amplitude as well as 
with Re. From Fig. 27 it is evident that for X/L = 0.20 and 
Re = 500 the back flow becomes stronger as well as more 
extensive. Moreover for higher X/L and Re, s^arated flow 
occurs in the coverging portion of the channel as well. 

From Figs. 14, IS, 22, and 26 we can make out the 
behaviour of the dimensionless fluid enthalpy for the same 
combinations of X /L and Re respectively. Prom these figures 
we notice that as the wall amplitude X/h is increased or the 
Reynolds number Re is decreased the difference b^ween the 
fluid enthalpy and the wall enthalpy decreases substantially*. 

Out of these two parameters the non-dimensional parameter 
X/L has the strongest effect in reducing this difference. 
Therefore, the length required to reach the wall enthalpy in 
the flow field decreases with Increase of X /L and with decrease 
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of Ro.* Prom a caroful insp_'Ction of samo f igur :s it. is 

soon that a point of infloction app-sars in tho onthalpv profil.-s 
in tho. Separated flow region whereas in the non-separatad 
region the enthalpy profile is almost parabolic- 

Fig- 23(a,b) and 29(a,b) show the dimensionless pressure 
distribution in the dOTrastream direction for the same X/L and 
Re- A simple analysis of ^igs» 23 and 29 mahes it possible 
to express this distribution as 

P(X) = -f(x) + p'(x) 

where p'(x) behaves in the periodic fashion from cycle to cycle 
and f(x) is nearly linear in the developing region and linear 
in the fully-developed region- The percycle pressure drop P 
in the fully developed region is same for each cycle- e«g- for 
X/L = 0,.10 and Re = 100- AP = -0.090. This per-cyclc pressure 
drop increases with increase of wall amplitude.* The variation 
of pressure in the Y-direction is negligible and has therefore 
not been shown* 

4.. 2 For Pip e 

Figs. 30 to 32, Figs- 34 to 36, Figs.- 33 to 40, Figs. 42 
to 44, and Figs* 46 to 43 show the behaviour of the dimensionless 
velocity components U and V and velocity vector U for X/L = 0.10 
and Re = 100, X/L = 0.10 and Re = 500, xA = 0-2° 

X/L = 0.20 and Re = 500, and x/L = 0.25 and Re = 100 respectively. 
Similar to that for the symmetric channel, the negative V- 
component at the minimum cross-section and the positive 
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V “Component at the maximiam cross-section increases with Re* 

In this case the effect of Re on the separated flow is stionger 
than that in the case of the symmetric channel* 

Here the flow begins to s^arate even at X/L = 0.*10 and Re = 100 
(of ■‘^ig* 30 ) since 9U/3Y -*■ 0 at some locations, e*g^ at X = 3*94 
near the wall* VJhen the Reynolds number increases from 100 to 
500 with the wall amplitude held constant, a slightly larger 
back flow than that in the case of the symmetric channel is 
observed at some locations, e*g* at X = 2*73, 3*94 in Fig* 34* 
Here also, as in the symmetric channel, a similar conclusion 
can be drawn about the effect of X/L and Re on separated flow 
by conparing the Pigs* 30, 34, 33, 42, and 46 * It is also 
observed that the back flow in this case is more stronger and 
extensive than that for the symmetric channel* 

Figs.* 33, 37, 41, 45, and 49 depict the dimensionless 
enthalpy H for the same ccmbinations of X/L and Rs.* Here also, 
the enthalpy is affected by the X/L and Re in a similar way 
as in the case of the symmetric channel* However, the length 
i-equired to reach the wall value is analler than that for the 
symmetric channel,. The point of inflection in the enthalpy 
profiles is also observed in the separated flow regions* 

Figs* 50 to 52 r^ resent the distribution of non- 
dimensional pressure P in the X-direction for the same values 
of X/L and Re. The pressure P(x) also varies in a similar 
fashion as in a symmetric channel* However the per-cycle 
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pressure drop A p in this case is higher for X /L := 0»10 

and Re = 100, AP = “0*243# As earlier, in the fully 
developed region is the same for each cycle and increases with 
increase of wall amplitude# Again th-^ variation of pressure 
in the Y-direction is negligible and has not been shown# 

It appears that the circular geometry enhances various 
effects in comparison to the p 1 an a-g geometry. 

4*3 For non“s^.mimo±ric Channel 

Figs# 53 t.o 55, Figs# 57 to 59, Figs# 61 to 63, and 
Figs# 6 5 to 67 show the behaviour of the dimensionless velocity 
components U and V and the velocity vector U for X/h = 0#10 
and Re = 100, X/L = 0#10 and Re = 500, X /L = 0,20 and Re = 100, 
and X/L = 0#20 and Re = 500 respectively# A comparison 
between Figs# 54 and 53 shows that the negative V—component 
at X = 15#50, 17 #50, etc# increases with increase of Reynolds 
number.# This is due to the effect of inertia [ 3 ]# From Figs# 
55, 57, 6 3, and 67 it is observed that the resultant velocity 
U is retarded near the lower boundary and accelerated near 
the upper boundary when the flow is in the downward direction 
but this behaviour is reversed vficn the flow is in the upward 
direction# Due to this behaviour flow separates near the 
upper boundary in one section whereas it s^aratos near the 
lower boundary in the adjacent section# It is also seen frm 
the same figures that the Reynolds number has a significant 
effect on the above phenomenon#- 
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Paxiin Pigs.# 53/ 57# and 61 it is observed that for 
X/L = 0#10 and Re = 100, 500, and for X/L =0.20 and Re = 100 
there is no s^aration in the flow field. But when Re increases 
to 500 with X/L = 0.20/ a anall back flow near the wall is 
observad. The effect of Re and X/L on the s^arated flow is 
similar to that in the symmetric channel. 

Pigs.. 56 / 60/ 64/ and 68 depict the non-dimensional 
enthalpy H for the same X/L and Re. Here alsO/ as in the 
s\mmetric channel or pipe/ the oathalpy is affected by X/L 
and Re in a similar manner. Point of inflection Is also observed 
in the separated region. 

Pigs. 69(a/b) and 70(a, b) show the distribution of 
non-dimensional pressure P in the X-direction for the same 
values of X/L and Re. Here# the pressure P varies almost in a 
linear fashion. Por X/L = 0.20 and Re = 500/ slight flactuation 
is observed in the distribution. These flactuation increases 
with increase of X/L and Re. Per-cycle pressure drop AP in 
the fully developcid region is same/ e#g. for X/L = 0.10/ and 
Re = 100, Ap = “0.283 and it increases with increase of 
amplitude. In the Y-direction a slight variation of pressure 
is observed in the developing region whereas in the fully- 
developed region this variation is negligible as shown in 
Pig. 71 for X/L = 0..20 and Re = 100. 
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PRESSURE DISTRIBUTION ALONG X FOR A SYN . CHANNEL 

WITH LA,NBDA=0.20 AND R^=I00 
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FIG. 51a PRESSURE DISTRIBUTION ALONG X FOR A PIPE 
WITH LANBDA*0 20 AND R«-^00 
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FIG. 52 PRESSURE DIS I RlBUTION ALONG X FOR A PI^^E 
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14.73 


16.17 


1 6. 83 


18 . 17 


15.50 


‘ 7.50 


FIG. 56 ENTHA.LPY PROFILES FOR A.NON-SYM. CH.4NNEi 

WITH L4MBDA=0 10 AND Re= ! 00 . 






[:vD=i .0 




I —i 

CVD=1 .0 



58 V-VELOCITY PROFILES FOR A NON-SYM . CHANNlL 
WITH LAMBDA=0 . 1 0 AND Re=50e 




k MMM ij\ A M AMAA 


134 


I — * 

ci^D-l 0 










138 



3 33 


C03=t . 










E.69 


7.39 


9. 19 


I 

CHD=1 . 



FIG. 68 enthalpy profiles FOR ^ NOK-SYKi . CHANNEL 

LITH L/\^13DA=e . 20 A^4D R^=500 



144 


4 



08 ^ 


rj.u# oya 

PRESSURE DISTRIBUTION ALONG X FOR A NON-SYM CHANNEL 

WITH LANBDA=-0 <0 AND R«-{00 



-2 48 -i 
-2 68 -^ 

-2 80 -^ 
at? J 


-o 48 A 
-2 68 -: 
-2 68 - 
-4 08 - 


FIG. 69b. 

PRESSURE DISTRIBUTION ALONG X FOR A nok-sYK 

WITH LAMBDA=0 "0 AND Re-500 


CHANNEL 



PRESSURE DISTRIBUTION ALONG X FOR A NON-SYM CHANNE 

WITH LAMBDA=0 20 AND Re=l00. 



PRESSURE DISTRIBUTION ALONG X FOR A NON'^YM CHANNEL 

WITH LAHBDA=0.20 AND Re=500 . 


CPD=1 .0 


0 50 ^ 


2 5l 


1 . 87 


3.33 


. I 0 4 88 



5.69 


7.39 


LPD=1 .0 


9.19 



ZP3=\ .0 


6.17 16.^ 


•'e. 17 


5.50 


17.50 


PRESSURE DISTRIBUTION ALONG Y FOR A NON-SYM . CHANNEL 

WITH LAMBDA=0.20 AND R^=100. 



Chapter 5 


CONCLUSION 

In the present ^^ 7 orh an analysis has been performed for 
the solution of two-dimensional laminar/ viscous flow through 
s;\mmetric as well as non-symmetric channels and through 
converging and diverging pipes. Suitable non-orthogonal 
algebraic coordinate transformations are used that map the 
physical irregular domain into a rectangular one. The 
conservation equations are derived on a control volume basis 
that ensures global as well as local conservation of mass, 
momentum and energy.# Also, velocity components U and U ^ are 
used as the primary dependent variables in the momontTam equations. 
Since those compon-onts are normal to the control volume faces, 
evaluation of the mass flow (P) and conv.ictive fluxes (d) 
passing through the faces is facilitated. 

In order to obtain the solution, the discretized 
conservation exqaations are solved using the SMPLEC algorithm 
[12] which reduces solution cost reasonably. Further reduction 
in the solution cost has been achieved by using the methodology 
sat forth in [12] for the solution of pressure correction 
equation. The computer code developed is extremely goaeral 
in that it can bo easily modified to consider turbulent flow 
and variable fluid properties. 

Results found for a fixed Pr and for various values 
of Re and ^/L show the following : 



In all the cases under consideration/ the separated 
flow region grows with increase of Re and X/L.. 

In the case of a symmetric channel or pipe/ the 
separated flow region extends mainly in the diverging 
portion but for higher Re and X/L it occurs in the 
converging portion as well. 

In the case of a non-symmetric channel/ the resultant 
velocity U is retarded near the lower boundary and 
accelerated near the upper boundary when the flow is 
in the downward direction but this behaviour is 
reversed when the flow is in the upward direction# 

Due to this flow si^arates near the upper boundary in 
one section while it s<:5)arates near the lower boundary 
in the adjacent section. 

In the case of a symmetric channel or pipo/ variation 
of pressure in the Y~dir action is negligible but in 
the case of a non-symmetric channel this variation is 
percv^tiblo in the developing region. 

The distribution of pressure in the X-dir action 
can be ej^ress as 

P(x) = -f(x) + p''( x) 

where p^(x) behaves in a periodic fashion from cycle 
to cycle/ and f(x) is nearly linear in the developing 
region and linear in the fully developed region for a 
symmetric channel or pipe. For a non-symmetric channel/ 
however/ P(X) is almost linear throughout. 



Per-cycle pressure drop Ap in the fully developed 
region is constant for each cycle and it increases with 
increase of Re and X/L^ 

A point of inflection is observed in the enthalpy 
profiles in the separated region whereas in the non- 
separated region this distribution is almost parabolic. 
In the case of symmetric channel or pipe, the negative 
V-cOTiponent at the minimum cross-section and positive 
V-component at the maximum cross-section increases with 
Re. This is due to the effect of inertia [3]. 

In general, the circular ga^metry enhances the various 
effects compared to the plane geometry. 



APPENDIX 


Discretization of Domain Using and F 2 Factor 

Here the factorsF^ and Fj are used for the purpose 
to generate a nonmniform grid soacing in both the directions.. 
VThere 

Pj as Expansion factor for grid stqj -length in the 
V -direction 

F 2 = Contraction or ejq^ansion factor for grid step- 

laigth in the g -direction, < 1 for symmetric 

channel or pipe and P 2 2^ 1 for non-Symmetric Cnannel. 


i) Discretization in the t)— direction 



{NXj-2) 

XU(NX^) Fj + ... + Pj 







XUCNXj) 



(XPIPE - 4 ) 



(NX--1) 



where f 


1 



ATJ^ VNXj. 

ii) Discretization in the | -direction 

a) For symmetric channel or pipe 

NyMi 
NyMl 

YV( 1 ) = 0 ,0 
YV( 2) = A5 ^ 

YVC3) = Agjj + 

3 

Z 
1 

W(N-adl) sAljj ^ '^(Nlii'12) 



as + Ag ^+( l+Pj+a. 


1 - F 


,M'yM 4 


*= 2 A§ jj + 


1-F. 


Let AS„ -•§ 


RPIPE , t At 


« f 2 A^i 


A$^ as RPIP^fj 


where 


pN:«iM 5 



# 


RPJPE 

yV(.NyMi.j 


All 
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b) For non-siTCmetric channel 


N/Ml 

NYMZ 


__ 


^^yD2-^^y t- 

2. 






♦1 ^ 


YV(NYD2) = (1 + 2+^2 ■*■•*•■** 

= 2(1+P2+^'2 +•'*+ p(NYD2-2)) 


DPIPE s 2 


p(NyD 2 .l) . 1 




A A|, = DPIPYf, 


(NTO2-1) . 1 
v^here £2 = 2 “—j; — 
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